

*******************************************************************
*********** Sovereign Bond Returns and Bond Prices - Overview  ****
*******************************************************************

/*

This do file produces the graphs and tables on sovereign bond returns and prices 
in Meyer, Reinhart and Trebesch "Sovereign Bonds since Waterloo" (QJE). 
The file is structured in 3 main parts:

I.) The codes in the preamble prepare the dataset for the subsequent analysis:
	1.1.: bond level analysis: restrict the sample, compute excess and yearly returns
	1.2.: preparatory work to construct country-portfolio returns 
	1.3.: preparatory work to construct global-portfolio returns
	
II.) Produces the main figures

III.) Produces the main tables

*/

*******************************************************************
****************** Part I: Preambel *******************************
*******************************************************************

*run the folowing codes if programs not installed 

* Standareturn_nomize settings accross users: 
* ieboilstart, version(15.0) 

* ssc install egenmore
* ssc install rangestat

* net install dm79, from("http://www.stata.com/stb/stb56") 
clear all
set more 1
capture log close


* Define format for output files 
	global format1 putexcel A1:Z200, fpattern(solid ,white  ,white) font(timesnewroman, 11) overwritefmt
	global format2 fpattern(solid ,white ,white) font(timesnewroman, 11)
		
* =======================================================
* 1.1 PREPARE Returns data  -  bond level  
* =======================================================
use returns_prices_1820_2016.dta, clear


* =======================================================
* Compute returns
* =======================================================
tsset bond_ID tm

* Calculate historical returns

* Adjust for default: adjust coupon payments due to defaults/partial payments
    gen coupon_spliced = coupon
	replace coupon_spliced = coupon * conversion_princ_synth if EMBIGsource == 0 & conversion_princ_synth ~=.
	gen coupon_adjust=(coupon_spliced*default_coupon_sharepaid)/12 if EMBIGsource == 0

* Adjust prices 
	gen price_spliced = price
	replace price_spliced = default_oldbond_spliced if default_oldbond_spliced ~= . & EMBIGsource == 0
	replace price_spliced = price * conversion_princ_synth  if conversion_princ_synth ~= . & EMBIGsource == 0
	 
* calculate monthly total nominal returns
	bys bond_ID (tm): replace return_nom=(price_spliced+coupon_adjust+default_cash)/L.price_spliced -1 if EMBIGsource == 0
	
* calculate monthly total real returns
	merge m:1 tm using "inflation.dta"
	drop if _merge==2
	drop _merge

	replace return_real = [((1 + return_nom) / (1 + inflation_tm_UK)) - 1] if currency == "GBP" & EMBIGsource == 0
	replace return_real = [((1 + return_nom) / (1 + inflation_tm_US)) - 1] if currency == "USD"  & EMBIGsource == 0
	drop inflation_tm_US inflation_tm_UK 
	
* Restrict EMBIG-Data 
	keep if (EMBIGmonthly_inclusion == "Yes" & EMBIGsource == 1)| EMBIGsource == 0
	
	
	save "returns_prices_1820_2016mod.dta", replace
	
* Restrict data 
	drop						 if inrange(y,1974,1994) // drop illiqubond_ID period
	replace return_real =. 		 if default_combined == 1 // drop double entries due to the linkage of prices in restructing events (old prices with new ones in the case of a bond exchange) 
	replace return_real =. 		 if default_bond_drop == 1 // drop double entries due to the linkage of prices in restructing events (old prices with new ones in the case of a bond exchange) 
	drop 						 if return_real == .  // drop missings
	drop 						 if amount_issUSD == . // There are seven bonds for which we were not able to find reliable data on amounts issued
	bys bond_ID y: egen countmonths = count(m)          // count month per year and Instrument for which we have returns data 
	drop 						 if countmonths <= 11 // restrict sample to full 12-month returns coverage per year
	drop countmonths
	
* Compute monthly excess returns
	merge m:1 tm using "other_asset_classes_returns.dta", ///	
	keepusing(tm y m returns_UKlongbond_real returns_USTreasury_real returns_USbills_real returns_UKbills_real)


	gen excessreturn_real 		= return_real -  returns_UKlongbond_real  if currency =="GBP"   // compute excess returns over US/UK 10year bonds
	replace excessreturn_real 	= return_real -  returns_USTreasury_real  if currency =="USD"
	gen excessbillsredeal 		= return_real -  returns_UKbills_real  	  if currency =="GBP"   // compute excess returns over US/UK 3- month treasury bills
	replace excessbillsredeal 	= return_real -  returns_USbills_real	  if currency =="USD"
 

* Compute yearly real bond returns accoreturn_noming to the equation 4 in the
* paper. 
* First, we generate natural logarithmic total by year and bond bond_ID. 
* Second, we create the returns 

	bys bond_ID y: egen lryreal = total( ln( return_real+1 ) ), missing  
	gen ryreal = exp( lryreal ) - 1 
	drop lryreal
	
* save dataset
	saveold loadData, replace


* =======================================================
* 1.2 PREPARE Returns data  -  country level 
* =======================================================

* Load data 
	use loadData, clear

* Compute monthly weighted country portfolio returns	
	#delimit ;	
	collapse 
		(mean) price return_nom return_real excessbillsredeal excessreturn_real 
		(max) default_nr_country default_serial 
		(lastnm) code
		[aw = amount_issUSD], by(country y m tm)
	;
	#delimit cr

* Restrict data 
	drop if return_real == .							// restrict sample to full 12-month returns coverage per year
	bys country y: egen countmonths = count(m)
	drop if countmonths != 12
	drop countmonths

* Compute yearly country portfolio returns
	foreach var in 	return_nom  return_real excessbillsredeal excessreturn_real{
	bys country y: egen lryreal = total( ln( `var'+1 ) ), missing
	gen `var'y = exp( lryreal ) - 1 
	drop lryreal	
	}	
	
* save country-level dataset
	saveold loadDatac, replace
	
	
* =======================================================
* 1.3 PREPARE Returns data  -  portfolio - time series level 
* =======================================================
* Load data 
	use loadData, clear

* Compute monthly global portfolio returns
	#delimit; 
	collapse (mean) return_nom return_real
		excessBonds = excessreturn_real 
		excessBills = excessbillsredeal 
		 [aw = amount_issUSD], by(y m tm)
	;
	#delimit cr 

* Compute yearly global portfolio returns, see Equation 4 in the paper
	
	bys y: egen lryreal = total( ln( return_real + 1 ) ), missing //real 
	gen ryreal = exp( lryreal ) - 1 
	drop lryreal

	bys y: egen lryreal = total( ln( return_nom + 1 ) ), missing // nominal 
	gen ry = exp( lryreal ) - 1 
	drop lryreal

	bys y: egen lryreal = total( ln( excessBonds + 1 ) ), missing // excess bonds
	gen ryrealExcess = exp( lryreal ) - 1 
	drop lryreal

* save time series portfolio returns data
	saveold loadDataTS, replace

* =======================================================
* 1.4 PREPARE other assets classes
* =======================================================
	use other_asset_classes_returns, clear
	
	foreach var of varlist returns*{
	local asset = substr("`var'", 9, .)
   	bys y: egen lryreal = total( ln( `var' + 1 ) ), missing 
	gen ry`asset' = exp( lryreal ) - 1 
	drop lryreal
	}

	collapse  ry* , by(y)
	
	gen rySafe_real = ryUKlongbond_real if y <1919
	replace rySafe_real = ryUSTreasury_real if y >=1919
	
	save other_asset_classes_returns_yearly, replace

* =======================================================
* 1.5 PREPARE default window
* =======================================================	
	use haircuts_final_sample, clear 
	
	keep if default_date ~=. 
	
	* Create default winodw 
	expand 72
	bys case_ID_final default_date: gen count = _n
	gen tm = default_date
	format %tm tm 
	replace tm = tm + count - 1
	drop count 

	* 2 years back
	expand 24, gen(expand)
	replace tm = default_date if expand == 1

	bys case_ID_final default_date tm expand: gen count = _n
	drop if count > 25 & expand == 1
	replace tm = tm-(count-1) if expand == 1 & count>=1
	drop if count ==  1 & expand == 1
	drop count expand 

	gen y=year(dofm(tm))
	gen m=month(dofm(tm))

	gen defaultSTART = 1 if tm == default_date
	gen defaultAROUND = 0 if defaultSTART== 1

	order country case_ID_final default_date tm defaultSTART defaultAROUND
	sort case_ID_final tm

	global lag = 72
	forvalues t = 1/$lag {
		bys case_ID_final (tm): replace defaultAROUND = `t' if defaultAROUND[_n-`t'] == 0
	}
	global lag = 24
	forvalues t = 1/$lag {
		bys case_ID_final (tm): replace defaultAROUND = -`t' if defaultAROUND[_n+`t'] == 0
	}	
	
	save default_window, replace
	
	
*******************************************************************
****************** Part II: FIGURES *************************************
*******************************************************************



* =================================================================
* Figure 3: Bond price sample: Coverage across time
* =================================================================
* Load returns data 
	use returns_prices_1820_2016, clear
		
* Count active bonds/countries	
	drop if price ==. 
	bys y: egen countries = nvals(country)
	bys y: egen bonds = nvals(bond_ID)

collapse bonds countries, by(y)

* tidy up 
	set obs `=_N+1'
	replace y = 1815 if y ==.
	tsset y 
	tsfill
	replace countries = 0 if countries == .
	replace bonds = 0 if bonds == .
	
* Moving average for number active countries/bonds
	tssmooth ma movC = countries , window(5 1 0)
	tssmooth ma movB = bonds , window(5 1 0)

* Produce and export figure 3
	#delimit; 
		twoway 
			(line movC y, cmissing(n) lwidth(medthick) color(gs7) lpattern(shortdash) yaxis(1)) 
			(area movB  y, cmissing(n)  color(olive) yaxis(2)) 
			,
			ylabel(0 (5) 60, nogrid angle(0)  axis(1) labsize(medium)) 
			ylabel(0 (50) 400, nogrid angle(0)  axis(2) labsize(medium)) 
			xlabel(1815 (20) 2015) 
			ytitle("Countries covered", axis(1) size(medlarge))  
			ytitle("Active bonds", axis(2) size(medlarge))  
			xtitle("") 
			graphregion(color(white)) 
			legend( symysize(*.5) symxsize(*.5) region(color(white)) col(1) 
			label(1 "Active bonds outstanding (right axis)") 
			label(2 "Countries covered (left axis)") size(medlarge) order(2 1))
	;
	#delimit cr 

	graph export "Figure3.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure3.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure3.tif", as(tif) replace
	
* =====================================================================
* Figure 4: Sovereign bonds in default, 1815-2016
* =====================================================================
* Load main dataset 
use returns_prices_1820_2016.dta, clear

* Count bonds
	bys tm default_bond: egen nvals = nvals(bond_ID)
	
	collapse nvals [aw = amount_issUSD ], by(default_bond tm y m)

	reshape wide nvals ,  i(tm y m)  j( default_bond)
	tsset tm
	tsfill
	replace m = month(dofm(tm))
	replace y = year(dofm(tm))


	replace nvals1 = 0 if nvals1 == . 
	replace nvals0 = 0 if nvals0 == . 


	gen share1 = nvals1/(nvals0+nvals1)*100
	gen min = 0 
	gen max = 100

	replace share1 = . if y>1980 & y<1994

* Produce and export figure 4 
	#delimit;
	twoway 
	(rarea min share1 y, bcolor(maroon) cmissing(n))
	(rarea share1 max y , bcolor(gs14)  cmissing(n)),
			graphregion(color(white)) xsize(10)
			ylabel(#5, nogrid angle(0) labsize(large) )  
			xlabel(1815 (20) 2015, labsize(large))
			xtitle("")   
			ytitle("Share in default, in %" "(weighted by principal amounts)", size(vlarge)) 
			legend( size(vlarge) region(color(white)) 
			rowgap(*0.0001) bmargin(tiny) symy(*0.5) symx(*0.5)
		label(1 "Bonds in default (full or partial), share weighted by principal amounts") 
		label(2 "Not in default, share weighted by principal amounts") 
			order(1 2 ) row(2))
	;
	#delimit cr

	graph export "Figure4.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure4.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure4.tif", as(tif) replace
	
* =====================================================================
* Figure 5: Distribution of bond returns and bond prices, 1815-2016
* =====================================================================
* Load return time series data
use loadDataTS, clear

collapse (mean) ryreal ry ,by(y)

replace ryreal = ryreal*100

* Compute and save locally statistics for figure 5 panel A
	sum ryreal ,det
	loc t1 : di %4.2f  r(skewness)
	loc t2 : di %4.2f  r(mean) 
	loc t3 : di %4.2f  r(p50) 
		

* Produce and export figure 5 panel A 	
	#delimit ;	
	histogram ryreal , normal percent  bin(30)
		normopts(lcolor(dknavy)) color(gs12) lwidth(0.1) 
		addplot(pci 0 `r(mean)' 15 `r(mean)', color(gs0*0.8))
		ylabel(0 (2.5) 15,nogrid angle(0) labsize(medlarge)  format(%4.0f))
		xlabel(, nogrid angle(0) labsize(medlarge) format(%4.0f)) 
		xtitle("") ytitle(, size(large))  
		graphregion(color(white))
		legend( size(large) symysize(*.3)  symxsize(*.3) region(color(white)) row(2) 
		label(1 "Total returns (yearly, real, arithmetic avg.)") 
		label(2 "normal distribution") 
		label(3 "") label(4 "") order(1 2))
		text(14 50 "Mean: `t2'" "Median: `t3'" "Scewness: `t1'" , lcolor(none) fcolor(none) 
		margin(small) box justification(left))
	 ;
	#delimit cr		
		
	graph export "Figure5_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure5_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure5_PanelA.tif", as(tif) replace

* ------------------------------------------------------------------------------

* Panel B: Load bond level returns data for figure 5 panel B
use loadData, clear

collapse ryreal = ryreal (lastnm) price (max) default_country, by(country bond_ID y)	

replace ryreal = ryreal*100

* Compute and save locally summary stats 
	sum ryreal if default_country~= 1 & ryreal <= 100,det
		loc t1 : di %4.2f  r(skewness)
		loc t2 : di %4.2f  r(mean) 
		loc t3 : di %4.2f  r(p50) 

	sum ryreal if  default_country== 1 & ryreal <= 100,det
		loc t6 : di %4.2f  r(skewness)
		loc t4 : di %4.2f  r(mean) 
		loc t5 : di %4.2f  r(p50) 

* Produce and export figure 5 panel B 
	#delimit;
	twoway 
	(hist ryreal if default_country~=1 & ryreal<=100, percent lcolor(gs11) lwidth(tiny) bin(30) fcolor(gs11)) 
	(hist ryreal if default_country==1 & ryreal<=100, percent fcolor(none) lcolor(cranberry) lwidth(tiny) bin(30)), 
	graphregion(color(white))
		title("") ytitle( ,size(large))
		xlabel(-100 (100) 100, angle(0) nogrid labsize(medlarge)) 
		ylabel(0 (5) 30, angle(0) nogrid labsize(medlarge) ) 
		xtitle("Total returns (yearly, real, arithmetic avg.)", size(large))
		legend(label(2 "Default years") label(1 "Non-default years") 
		order(1 2) region(color(white)) symy(*1.8) symx(*0.5) col(2) )
		text(25 75 "{bf:Non-default years:}" "Mean: `t2'" "Median: `t3'" "Skewness: `t1'" " " "{bf:Default years:}" "Mean: `t4'" "Median: `t5'" "Skewness: `t6'", lcolor(none) fcolor(none) 
		margin(small) box justification(left))
	;
	#delimit cr

	graph export "Figure5_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure5_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure5_PanelB.tif", as(tif) replace

* ------------------------------------------------------------------------------
* Panel C

* Compute and save locally summary stats for figure 5 panel C 
	sum price if default_country ~= 1 & price <= 150,det
		loc t1 : di %4.2f  r(skewness)
		loc t2 : di %4.2f  r(mean) 
		loc t3 : di %4.2f  r(p50) 

	sum price if  default_country == 1 & price <= 150,det
		loc t6 : di %4.2f  r(skewness)
		loc t4 : di %4.2f  r(mean) 
		loc t5 : di %4.2f  r(p50) 

* Produce and export figure 5 panel C  
	#delimit;
	twoway 
	(hist price if  default_country ~= 1 & price<=150 ,  percent  lcolor(gs11) lwidth(tiny) bin(30) fcolor(gs11)) 
	(hist price if  default_country == 1  & price<=150 ,  percent fcolor(none) lcolor(cranberry) bin(30) lwidth(tiny)), 
	graphregion(color(white))
		title("") ytitle(, size(large))
		ylabel( ,angle(0) nogrid labsize(medlarge)) 
		xlabel(, nogrid angle(0) labsize(medlarge) format(%4.0f)) 
		xtitle("Bond prices (end of year)" , size(large))
		legend(label(1 "Non-default years") label(2 "Default years") order(1 2) 
		region(color(white)) symy(*1.8) symx(*0.5) col(2) )
		text(16 130 "{bf:Non-default years:}" "Mean: `t2'" "Median: `t3'" "Skewness: `t1'" 
		" " "{bf:Default years:}" "Mean: `t4'" "Median: `t5'" "Skewness: `t6'", 
		lcolor(none) fcolor(none) 
		margin(small) box justification(left) )
	;
	#delimit cr


	graph export "Figure5_PanelC.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure5_PanelC.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure5_PanelC.tif", as(tif) replace

* =====================================================================
* Figure 6: The bottom bin: Share of bonds with negative returns
* =====================================================================
* Load bond level returns data
	use loadData, clear
	tsset bond_ID tm

* Prepare and tbond_IDy up dataset for figure 6	
	forvalues i = 1/10{
	sort bond_ID tm 
	local ii = ((12*`i')-1)
	dis `ii'
	gen double lret1 = log(return_real+1)
	rangestat (sum) lret1 (count) lret1, interval(tm -`ii' 0) by(bond_ID)
	gen wanted`i' = exp(lret1_sum) - 1 if lret1_count == `ii'+1
	drop lret1 lret1_sum lret1_count
	gen majorloss`i' = 1 if wanted`i'<= -0.3&wanted`i' ~= .
	replace majorloss`i' = 0 if wanted`i'>-0.3&wanted`i' ~= .
	gen hugeloss`i'=1 if wanted`i'<= -0.6&wanted`i' ~= .
	replace hugeloss`i' = 0 if wanted`i'>-0.6&wanted`i' ~= .
	gen wantedd`i' = 1 if wanted`i'<0 & wanted`i'~=.
	replace wantedd`i' = 0 if wanted`i' >=0 &wanted`i' ~= .
	replace wanted`i' = wanted`i'*100
	}

	preserve
	tempfile temp1
	tempfile temp2
	clear 
	gen n=1
	save `temp1', replace
	save `temp2', replace
	restore


	foreach stats in sum count{
	preserve
	collapse (`stats') wantedd* majorloss* hugeloss*
	gen code = 1
	reshape long wantedd majorloss hugeloss, i(code) j(n)
	drop code
	rename wantedd `stats'
	rename majorloss loss_`stats'
	rename hugeloss hugeloss_`stats'
	merge 1:1 n using `temp1', update
	drop _merge
	save `temp1', replace
	restore 
	}


	foreach stats in sum count{
	preserve
	drop if y<1990
	collapse (`stats') wantedd* majorloss* hugeloss*
	gen code=1
	reshape long wantedd majorloss hugeloss, i(code) j(n)
	drop code
	rename wantedd `stats'
	rename majorloss loss_`stats'
	rename hugeloss hugeloss_`stats'
	merge 1:1 n using `temp2', update
	drop _merge
	save `temp2', replace
	restore 
	}



	clear 
	use `temp1' 
	gen code = 2

	append using  `temp2' 

	gen fre1 = (sum/count)*100
	gen freloss1 = (loss_sum/count)*100
	gen huge1 = (hugeloss_sum/hugeloss_count)*100

	replace code = 3 if code ==.

* ------------------------------------------------------------------------------
* Panel A: Produce and export figure 6 panel A

	#delimit;
	twoway 
		(line fre1 n if code == 2, color(gs4) lpattern(longdash)) 
		(line freloss1 n if code == 2, color(gs4)) 
		(line huge1 n if code == 2, color(cranberry) 
					ylabel(0 (10) 30, labsize(large) angle(0) nogrid format(%9.0f)) 
					xlabel(1(1)10, labsize(large)) 
					lwidth(medthick) 
					fintensity(200) 
					) 
			, title("", color(gs0)) 
			ytitle("Share of bonds in %", size(large)) 
			xtitle("Holding period in years", size(large)) 
			legend(rows(3) size(large) 
			label(1 "Share of bonds with returns < 0%") 
			label(2 "Share of bonds with returns < -30%") 
			label(3 "Share of bonds with returns < -60%") 
			region(color(white))   
			 symysize(5) symxsize(5)  rowgap(-3) linegap(-1) margin(zero) bmargin(tiny))  
			graphregion(color(white))  
	;
	#delimit cr


	graph export "Figure6_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure6_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure6_PanelA.tif", as(tif) replace

		
* ------------------------------------------------------------------------------
* Panel B: Produce and export figure 6 panel B

	#delimit;
	twoway 
		(line fre1 n if code == 3, color(gs4) lpattern(longdash)) 
		(line freloss1 n if code == 3, color(gs4)) 
		(line huge1 n if code == 3, color(cranberry) 
			ylabel(0 (10) 30,labsize(large) angle(0) nogrid format(%9.0f)) 
			xlabel(1(1)10, labsize(large) ) 
			lwidth(medthick) 
			fintensity(200) 
			) 
			, title("", color(gs0)) 
			ytitle("Share of bonds in %", size(large)) 
			xtitle("Holding period in years" , size(large)) 
			legend(rows(3)  size(large) 
			label(1 "Share of bonds with returns < 0%") 
			label(2 "Share of bonds with returns < -30%") 
			label(3 "Share of bonds with returns < -60%") 
			region(color(white)) rowgap(-3) linegap(-1) margin(zero) bmargin(zero)  
			 symysize(5) symxsize(5)) 
			graphregion(color(white))
	;
	#delimit cr


	graph export "Figure6_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure6_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure6_PanelB.tif", as(tif) replace


* =====================================================================
* Figure 7
* =====================================================================

* Load main dataset
	use returns_prices_1820_2016mod.dta, clear
	
* Restrict sample 	
	drop if inrange(y,1974,1994)
	replace return_real = . 			 if default_combined == 1
	drop if return_real == . 
	drop if amount_issUSD == .
	bys bond_ID y: egen countmonths = count(m)
	drop if countmonths <= 11 & y>1860
	drop countmonths

collapse (mean) return_real [w = amount_issUSD],by(default_serial y m tm)

* Compute returns by default category
	bys default_serial y: egen lryreal = total( ln( return_real+1 ) ), missing
	gen ryreal = exp( lryreal ) - 1 
	drop lryreal

* Add UK/US returns
	append using "other_asset_classes_returns_yearly.dta", keep(y rySafe_real)
	replace default_serial = 2 if default_serial == . 
	replace ryreal = rySafe_real if default_serial == 2


collapse (mean) ryreal, by( default_serial y)


* Define events
	gen events = 1 if inrange(y,1914,1918)
	replace events = 1 if inrange(y,1939,1945)
	
* Compute moving average	
	tsset default_serial y
	tsfill
	tssmooth ma mov = ryreal, window(5 1 5)
	keep if y >= 1815
	replace mov = . if inrange(y,1974,1994)
	replace mov = mov*100

* tbond_IDy up	
	quietly sum mov
	gen max = r(max) if events ~= .
	gen min = -10 if events ~= .
	
* Figure 
* ------------------------------------------------------------------------------
* Panel A: Produce and export figure 7 panel A
		
	#delimit;	
	twoway (area max y , cmissing(no) color(gs15)) 
			(area min y , cmissing(no) color(gs15)) 
			(function y = 0, range(1815 2015) clstyle(dash) clpattern(longdash)) 
			(line mov y if default_serial == 2 , cmissing(n) lwidth(medthick) 
				lcolor(gs7) lpattern(longdash)) 
			(line mov  y if default_serial == 1, cmissing(n) lwidth(medthick) 
			lcolor(cranberry)) , 
			ylabel(-10(5)15,nogrid angle(0) labsize(huge)) 
			xlabel(1815 (20) 2015, labsize(huge)) 
			ytitle("") xsize(12) 
			xtitle("") 
			graphregion(color(white)) 
			legend( symysize(*.5) size(huge) symxsize(*.5) 
			region(color(white)) row(1) 
			label(4 "UK/US 10y bonds") 
			label(5 "Serial defaulters") order(4 5))
			
	;
	#delimit cr
	graph export "Figure7_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure7_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure7_PanelA.tif", as(tif) replace
	
* ------------------------------------------------------------------------------
* Panel B: Produce and export figure 7 panel B
	
	#delimit;	
	twoway (area max y , cmissing(no) color(gs15)) 
			(area min y , cmissing(no) color(gs15)) 
			(function y = 0, range(1815 2015) clstyle(dash) clpattern(longdash)) 
			(line mov y if default_serial == 2 , cmissing(n) lwidth(medthick) 
				lcolor(gs7) lpattern(longdash)) 
			(line mov  y if default_serial == 0, cmissing(n) lwidth(medthick) 
			lcolor(gs0)) , 
			ylabel(-10(5)15,nogrid angle(0) labsize(huge)) 
			xlabel(1815 (20) 2015, labsize(huge)) 
			ytitle("") xsize(12) 
			xtitle("") 
			graphregion(color(white)) 
			legend( symysize(*.5) size(huge) symxsize(*.5) 
			region(color(white)) row(1) 
			label(4 "UK/US 10y bonds") 
			label(5 "All others") order(4 5))
	;
	#delimit cr
	graph export "Figure7_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure7_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure7_PanelB.tif", as(tif) replace


* =====================================================================
* Figure 8: Sovereign risk-return profiles: Returns and standard deviations, 1815-2016
* =====================================================================
* Load bond level returns data
	use loadData, clear

* Drop outlier 
	drop if country == "Romania" & y == 1952
	drop if country == "Romania" & y == 1959

collapse (mean) return_real (max) default_nr_country default_serial, by(country code y m tm)

* Compute yearly real country portfolio returns 
	bys country y: egen lryreal = total( ln( return_real+1 ) ), missing
	gen ryreal = exp( lryreal ) - 1 
	drop lryreal

collapse (mean) ryreal (max) default_nr_country default_serial, by(country code y)

* Add alternative asset classes
	append using "other_asset_classes_returns_yearly.dta", keep(y ryUSTreasury_real)
	replace code = "US bonds" 				if code == ""
	replace ryreal = ryUSTreasury_real 		if code == "US bonds"
	replace default_nr_country = 0 					if code == "US bonds"
	replace default_serial = 2 					if code == "US bonds" 

	append using "other_asset_classes_returns_yearly.dta", keep(y ryUKlongbond_real)
	replace code = "UK bonds" 			if code == ""	
	replace ryreal = ryUKlongbond_real		if code == "UK bonds"
	replace default_nr_country = 0 				if code == "UK bonds"
	replace default_serial = 2 				if code == "UK bonds" 

	replace country = code if country == ""


collapse (mean) ryreal (max) default_nr_country default_serial (sd) SDr= ryreal ///
		(count)  y=ryreal, by(country code)

* Tidy up 
	foreach var of varlist ryreal SDr {
	replace `var' = `var' * 100
	}

* Drop coverage with less than 10 years
	drop if y<10
	
* Compute p-value and t- statistics	for figure 8
	reg ryreal SDr 
	loc t : di %4.2f  _b[SDr]/_se[SDr]  // get the t-statistic
	loc p : di %4.2f 2*ttail(e(df_r),abs(`t'))   //Calculate p' value

* Produce and export figure 8 
		
	#delimit;
	twoway 
			(lfit ryreal SDr, lcolor(gs14))
			(scatter ryreal SDr if default_serial == 0  , mcolor(gs0) msymbol(smcircle) 
					mlabel(code) mlabcolor(gs0)) 
			(scatter ryreal SDr if default_serial == 1 , mcolor(cranberry) msymbol(smplus) 
					mlabel(code) mlabcolor(cranberry)) 
			(scatter ryreal SDr if code=="UK bonds", mcolor(blue) msymbol(circle) 
					mlabel(code) mlabcolor(blue) mlabposition(4) msize(medium) mlabsize(medsmall) mlabgap(1) ) 
			(scatter ryreal SDr if code=="US bonds", mcolor(blue) msymbol(circle) 
					mlabel(code) mlabcolor(blue) mlabposition(12) msize(medium) mlabsize(medsmall) mlabgap(2) ) ,		
			ylabel(,nogrid angle(0) labsize(medlarge)) 
			xlabel(0 (10) 50, labsize(medlarge)) 
			ytitle("Mean annual real return in %", size(medlarge)) 
			xtitle("Standard deviation", size(medlarge)) 
			graphregion(color(white)) 
			plotregion(margin(medlarge)) 
			legend( rows(1) symysize(*5) size(medlarge) symxsize(*5) label(3 "Serial defaulters")
			label(2 "All others") order(3 2)  bmargin(tiny) 
			region(color(white))) 
			text(0.2 3.2 "t-statistic:`t'" "p-value:`p'" , lcolor(gs0) fcolor(none) 
			margin(small) box justification(left))
			
	;
	#delimit cr		
	graph export "Figure8_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure8_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure8_PanelA.tif", as(tif) replace


* =====================================================================
* Figure 9: The role of credit history: Returns and default frequency, by country, 1815-2016
* =====================================================================

* Compute p-value and t- statistics for figure 9
	reg ryreal default_nr_country 
	loc t : di %4.2f  _b[default_nr_country]/_se[default_nr_country]  // get the t-statistic
	loc p : di %4.2f 2*ttail(e(df_r),abs(`t'))   //Calculate p' value


* Produce and export figure 9
	#delimit;
	 twoway 
			(lfit ryreal default_nr_country, lcolor(gs14))
			(scatter ryreal default_nr_country if default_serial==0 , mcolor(gs0) msymbol(smcircle) 
					mlabel(code) mlabcolor(gs0)) 
			(scatter ryreal default_nr_country if default_serial==1 & code~="TUR", mcolor(cranberry) msymbol(smplus) 
					mlabel(code) mlabcolor(cranberry)) 	
			(scatter ryreal default_nr_country if default_serial==1 & code=="TUR", mcolor(cranberry) msymbol(smplus) 
					mlabel(code) mlabcolor(cranberry) mlabposition(6)) 
			(scatter ryreal default_nr_country if code=="UK bonds", mcolor(blue) msymbol(circle) 
					mlabel(code) mlabcolor(blue) mlabposition(5) msize(medium) mlabsize(medsmall) mlabgap(3) ) 
			(scatter ryreal default_nr_country if code=="US bonds", mcolor(blue) msymbol(circle) 
					mlabel(code) mlabcolor(blue) mlabposition(5) msize(medium) mlabsize(medsmall) mlabgap(-1) ) , 
			ylabel(,nogrid angle(0) labsize(medlarge)) 
			xlabel(0 (1) 11, labsize(medlarge)) 
			ytitle("Mean annual real return in %", size(medlarge)) 
			xtitle("Number of external defaults", size(medlarge)) 
			graphregion(color(white)) 
			legend( rows(1) symysize(*5) size(medlarge) symxsize(*5) label(3 "Serial defaulters") 
			label(2 "All others") order(3 2)  bmargin(tiny) 
			region(color(white)))		
			text(1 10 "t-statistic:`t'" "p-value:`p'" , lcolor(gs0) fcolor(none) 
			margin(small) box justification(left))
			
	;
	#delimit cr		
	graph export "Figure8_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure8_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure8_PanelB.tif", as(tif) replace

* =====================================================================
* Figure 10: Bid-ask spreads and total returns, 1870-2016
* =====================================================================
* Load bond level returns data
	use returns_prices_1820_2016mod, clear
	
* Compute relative bbond_ID-ask prices
	gen BIDASKspreadbond_IDrel = (pask-pbid)/(0.5*(pask+pbid))*100 if pbid ~= . & pask ~= . & pask>=pbid

collapse BIDASKspreadCoutnriesREL = BIDASKspreadbond_IDrel return_real, by(tm y m)

* Tidy up 			
	tsset tm
	tsfill
	order y m tm  BIDASKspreadCoutnriesREL
	replace y = year(dofm(tm))
	replace m = month(dofm(tm))
	sort tm
	replace return_real = return_real*100 

* Create event dummy  	
	gen events = ""
	 replace events = "London Panic"             if inrange(y,1826,1830)
	 replace events = "Baring Crises"            if y == 1890
	 replace events = "Founders crash"           if y == 1880
	 replace events = "WWI" 			         if inrange(y,1915,1918)
	 replace events = "WWII"                     if inrange(y,1939,1945)
	 replace events = "Great Drepression"        if inrange(y,1929,1934)
	 replace events = "Global Debt Crisis 1980s" if inrange(y,1980,1987)
	 replace events = "Financial crisis"         if inrange(y,2007,2008)
	 replace events = "Founders crash"           if y == 1880
	 replace events = "Panic of 1907"            if y == 1907
	 replace events = "Wiener Börsenkrach"       if y == 1873
	 replace events = "defaults"                 if y == 1875 & m >= 10
	 replace events = "defaults"                 if y == 1876 & m <= 5
	 replace events = "defaults"                 if y == 1893 & m <= 7
	 replace events = "defaults"                 if y == 1891
	 replace events = "defaults"                 if y == 1892
	 replace events = "defaults"                 if y == 1920
	 replace events = "Russia default"           if y == 1998 & m >= 7
	 replace events = "Argentina default"        if y == 2001 & m >= 11
	 replace events = "Lehmann"                  if y == 2008 & m >= 9
	 replace events = "Indonesia"                if y == 1997 & m >= 7
	 replace events = "Greek restructuring"      if y == 2011 & m >= 10
	 replace events = "" 						 if inrange(y,1980,1990)
 
* Compute max/min for figure 9 panel A 
	gen max = 10 if events ~= ""
	gen min = -5 if events ~= ""

* Compute moving average for figure 9 panel A
	tssmooth ma sharp1movLIQ = BIDASKspreadCoutnriesREL, window(6 1 6)
	tssmooth ma sharp1movEXCESS = return_real, window(6 1 6)
		
* Figure 
* ------------------------------------------------------------------------------
* Produce and export figure 9 panel A
	#delimit;
	twoway 
		(area max tm if y < 1970 & y >= 1870 , cmissing(no) color(gs14)) 
		(area min tm if y < 1970 & y >= 1870, cmissing(no) color(gs14))
		(area sharp1movEXCESS tm if y < 1970 & y >= 1870, color(maroon*0.6)cmissing(no) lwidth(none))
		(line sharp1movLIQ tm if y < 1970 & y >= 1870, color(dknavy) cmissing(no) lwidth(medthick))
		, 
		xsize(12) graphregion(color(white) )
		xtitle("")
		xlabel(#20, format("%tmCCYY") labsize(vlarge))
		ylabel(-5 (5) 10, angle(0) nogrid labsize(vlarge)) 
		ytitle("" ,size(vlarge))
		legend( size(large) region(color(white)) symy(*1.8) symx(*0.5) row(1) 
		label(3 "Total real returns (monthly moving average, in %)") 
		label(4 "Bid-ask spreads (monthly moving average, in %)") 
		order(3 4)
		margin(medsmall) bmargin(medsmall))
		text(9 -632  "{bf:Panic}" "{bf:of 1907}", size(small))
		text(9 -1041  "{bf:Vienna}" "{bf:crash}", size(small))
		text(9 -526  "{bf:World}" "{bf:War I}", size(small))
		text(9 -200  "{bf:World}" "{bf:War II}", size(small))
		text(9 -337  "{bf:Great}" "{bf:Drepression}", size(small))
		text(9 -957  "{bf:Panic}" "{bf:of 1873}", size(small))
		text(-4 -1004  "{bf:Egypt}" "{bf:and}" "{bf:Ottoman Empire}" "{bf:default}", size(small))
		text(8.5 -820  "{bf:Baring}" "{bf:crisis}" "{bf:and}" "{bf:default}" "{bf:Argentina}", size(small))
		text(8.5  -470  "{bf:Russian}" "{bf:default,}" "{bf:Empire}" "{bf:break-ups}", size(small))
	;
	#delimit cr 

	graph export "Figure9_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure9_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure9_PanelA.tif", as(tif) replace

* -------------------------------------------------------------------------------------

* Compute max/min for figure 9 panel B
	capture drop max min
	gen max = 5 if events ~= ""
	gen min = -5 if events ~= ""

* Produce and export figure 9 panel B 	
	#delimit;
	twoway 
		(area max tm if y>1994 , cmissing(no) color(gs14)) 
		(area min tm if y>1994, cmissing(no) color(gs14))
		(area sharp1movEXCESS tm if y>1994, color(maroon*0.6)cmissing(no) lwidth(none))
		(line BIDASKspreadCoutnriesREL tm if y>1994, color(dknavy) cmissing(no) lwidth(medthick)), 
		xsize(12) graphregion(color(white) )
		xtitle("")
		xlabel(#20, format("%tmCCYY") labsize(vlarge))
		xtitle("")
		ylabel( -5  0 5, angle(0) nogrid labsize(vlarge)) 
		ytitle("")
		legend( size(vlarge) region(color(white)) symy(*1.8) symx(*0.5) row(1) 
		label(3 "Total real returns (monthly moving average, in %)") 
		label(4 "Bid-ask spreads (monthly moving average, in %)") 
		order(3 4)
		margin(medsmall) bmargin(medsmall) )
		text(4.8 465 "{bf:Russia}" "{bf:default}" "{bf:and}" "{bf:LTCM crisis}", size(small))
		text(5 500 "{bf:Argentina defaults,}" "{bf:Brazil crisis}", size(small))
		text(5 585 "{bf:Lehmann}" "{bf:bankruptcy}", size(small))
		text(5 452 "{bf:Asian}" "{bf: crisis}", size(small))
		text(5 620 "{bf:Greek}" "{bf:restructuring}" "{bf:announcement}", size(small))
	;
	#delimit cr 
	drop min max

	graph export "Figure9_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure9_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure9_PanelB.tif", as(tif) replace

* =====================================================================
* Figure 11: Bond returns around sovereign defaults
* =====================================================================
* Load default panel data set 
	use default_window, clear 
	
* Merge country portfolio returns	
	merge m:1 country tm y m using loadDatac, keep(match)
	drop _merge
	
* Select only deals with sufficient data 
	bys case_ID_final (defaultAROUND): egen count = count(return_real)
	gen keep = .
	replace keep = 1 if count == 96
	foreach ID in 2 48 51 123 198 219 251 60 62 136 189 198 199 287{
	replace keep = 1 if case_ID_final == `ID'
	}
	keep if keep == 1
	drop keep	

* Compute cumulative return	
	bys country case_ID_final (defaultAROUND): gen cumuReturn = log(1) if defaultAROUND == -24 
	bys country case_ID_final (defaultAROUND): replace cumuReturn = cumuReturn[_n-1]+(log(1+return_real)) if defaultAROUND>-24 
	replace cumuReturn = exp(cumuReturn)
	egen countR = count(cumuReturn), by(case_ID_final)


#delimit ;	
collapse 
(mean) cumuReturn
(p75)  cumuReturn75 = cumuReturn
(p25)  cumuReturn25 = cumuReturn
, by(defaultAROUND)
;
#delimit cr	
	  

* Tidy up 
	foreach var of varlist cumuReturn25 cumuReturn cumuReturn75{
	replace `var' = `var'* 100
	}

* Produce and export figure 10  
	#delimit;
	twoway	 
			(line  cumuReturn25 defaultAROUND , lwidth(medthick) lpattern(dash) color(gs8))
			(line  cumuReturn defaultAROUND , lwidth(medthick) color(gs0)) 
			(line  cumuReturn75 defaultAROUND , lwidth(medthick) lpattern(dash) color(gs8))
			(pcarrowi 165 50 150 55, mcolor(gs0) lcolor(gs0))
			(pcarrowi  75 50 60 55, mcolor(gs0) lcolor(gs0))
			(function y=100 , range(-24 80) color(gs10) 
			lwidth(medium) clpat(dash)) 
			, ylabel(50 75 100 125 150 175 ,nogrid angle(0)  format(%4.0f) labsize(medlarge)) 
			graphregion(color(white)) 
			ytitle("Cumulative total return index", size(medlarge)) 
			xtitle("Years around default", size(medlarge)) 
			text(117 68 "Average" "cumulative" "return", size(medium) place(e))
			text(170 40 "Upper 75%", size(medium) place(e))
			text(80 40 "Lower 25%", size(medium) place(e))
			xlabel(-24 "-2" -12 "-1" 0 "t" 12 "1" 24 "2" 36 "3" 48 "4" 60 "5" 72 "6" , labsize(medlarge)) 
			legend( rows(3)  size(medlarge) rowgap(*0.0001) bmargin(tiny) 
			symysize(*.5)  symxsize(*.5) region(color(white)) 
			label(1 "Start of holding period: 2 years before default")  
			order(1))
			
	;
	#delimit cr	

	graph export "Figure10.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure10.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure10.tif", as(tif) replace

* =====================================================================
* Figure 12: Returns around default: High vs. low haircut cases       
* =====================================================================
* Load default panel data set 
	use default_window, clear 
	
* Merge country portfolio returns	
	merge m:1 country tm y m using loadDatac, keep(match)
	
* Select only deals with sufficient data 
	bys case_ID_final (defaultAROUND): egen count = count(return_real)
	gen keep = .
	replace keep = 1 if count == 96
	foreach ID in 2 48 51 123 198 219 251 60 62 136 189 198 199 287{
	replace keep = 1 if case_ID_final == `ID'
	}
	keep if keep == 1
	drop keep	
	
* Get median haircut
	preserve
	bys case_ID_final: keep if _n == 1
	sum haircut_pv, det
	local med = r(p50)
	restore
		
* Compute cumulative return	
	bys country case_ID_final (defaultAROUND): gen cumuReturn=log(1) if defaultAROUND == -24 
	bys country case_ID_final (defaultAROUND): replace cumuReturn=cumuReturn[_n-1]+(log(1+return_real)) if defaultAROUND>-24 
	replace cumuReturn = exp(cumuReturn)

* Split sample by median haircut 	
	gen dumhigh = 1 if haircut_pv > `med' & haircut_pv~=.
	replace dumhigh = 0 if haircut_pv <= `med' & haircut_pv~=.
	drop if dumhigh ==.

	
#delimit ;	
collapse 
(mean) cumuReturn
(p75) cumuReturn75 = cumuReturn
(p25) cumuReturn25 = cumuReturn
, by(dumhigh defaultAROUND)
;
#delimit cr	

replace cumuReturn = cumuReturn*100


* Produce and export figure 11  
	#delimit;
	twoway	 
			(line  cumuReturn defaultAROUND if dumhigh == 1, lwidth(medthick)  color(gs0))
			(line  cumuReturn defaultAROUND  if dumhigh == 0, lwidth(medthick) color(gs9)) 
			(pcarrowi 120 25 110 35, mcolor(gs9) color(gs9))
			(pcarrowi 60 55 70 50, mcolor(gs0) lcolor(gs0))
			(function y=100 , range(-24 80) color(gs10) 
			lwidth(medium) clpat(dash)) 
			, ylabel(50 75 100 125 150 175 ,nogrid angle(0)  format(%4.0f) labsize(medlarge)) 
			graphregion(color(white)) 
			ytitle("Cumulative total return index", size(medlarge)) 
			xtitle("Years around default", size(medlarge)) 
			text(125 0 "Low/moderate" "haircut defaults" "(haircut <=42%)", size(medium) place(e))
			text(60 55 "High" "haircut defaults" "(haircut >42%)", size(medium) place(e))
			xlabel(-24 "-2" -12 "-1" 0 "t" 12 "1" 24 "2" 36 "3" 48 "4" 60 "5" 72 "6" , labsize(medlarge)) 
			legend( rows(3)  size(medlarge) rowgap(*0.0001) bmargin(tiny) 
			symysize(*.5)  symxsize(*.5) region(color(white)) 
			label(1 "High haircut defaults: total return index (year t-2=100)") 
			label(2 "Low/moderate haircut defaults: total return index (year t-2=100)")
			order(2 1))
			
	;
	#delimit cr	
	
	graph export "Figure11.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure11.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure11.tif", as(tif) replace



* =====================================================================
* Figure 13: Returns vs. haircuts across years - historical and modern sample
* =====================================================================
* Load haircut data
	use haircuts_full_sample.dta, clear

	gen defaulty = year(dofm(default_date))
	
* Keep only bond exchange
	preserve
	drop if bond_restructuring == 0 & year_agreement > 1973
	collapse haircut_pv [aw = debt_restructured_real], by(country defaulty)
	rename haircut_pv hc_defy
	gen ID = "7777" + string(_n)
	destring ID, replace
	tempfile temp
	save "`temp'", replace	
	restore

* Load bond level returns data 
	use loadData, clear
	drop _merge 

* merge with default data 
* Merge with default data. Note: During the merging process over 205,000 
* observations will be not matched because only a fraction of total sovereign
* bonds are restructured, resulting in some haircut. Therefore, haircut
* entries are the reason for this high level of not-matched data. Instead, 
* those over 205,000 not-matched entries indicate that during a given year in a
* given country, no haircut was inccured. Thus, we do not drop those not-
* matched entries and replace them with 0 whenever no entry for haircuts
* was given. Thus, we can make a sound aggregation of average annual haircut
* and returns data. 

	rename y defaulty
	merge m:1 country defaulty using "`temp'", update
	drop _merge
	rename defaulty y
	replace hc_defy = 0 if hc_defy == .
	
* Compute average weighted haircut per year 
	collapse (mean)  hc_defy [aw = amount_issUSD],by(y)	
	
	tempfile temp
	save "`temp'", replace

* Load time series returns data 
	use loadDataTS, clear
	
* Compute yearly portfolio returns 
	collapse ryreal=ryreal, by(y)

* Merge with yearly weighted realized haircuts 
	merge 1:1 y using `temp'

* Tidy up 
	replace  ryreal = ryreal*100
	gen d="Modern bond period" if inrange(y,1995,2016)
	replace d="Historical bond period" if inrange(y,1815,1973)

* Produce and export figure 12 
	#delimit ;	
		graph bar ryreal  
		hc_defy, over(d,
		relabel(1 `""Historical bond period" "(1815-1973)""'
		2 `""Modern bond period" "(1995-2016)""')
		gap(10) label(angle(0) labsize(large))) 
		ylabel(,nogrid angle(0) labsize(large))
		ytitle("Return/haircut in %", size(large)) 
		bar(1,color(gs8)) bar(2,color(maroon)) intensity(65) lintensity(10) 
		legend( rowgap(*0.01) colgap(*0.01)
		label(1 "Returns (real, avg. annual)") 
		label(2 "Losses (avg. annual haircut)") size(large) order(1 2) ) 
		graphregion(color(white)) 
		bargap(-20) 
		legend( symysize(*.5)  symxsize(*.5) size(large) 
		rowgap(*0.01) colgap(*0.01)
		region(color(white)) row(3)) 
	;
	#delimit cr

	graph export "Figure12.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure12.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure12.tif", as(tif) replace

* =====================================================================
* Figure 14:  Asset classes across 200 years: Risk and return
* Done together with Table 7 (see below the code following Table 7)
* =====================================================================


* =====================================================================
* Figure 15
* =====================================================================
* Load main dataset 
	use returns_prices_1820_2016mod.dta, clear
	
	
* Restrict
	drop if y>2016
	drop if inrange(y,1974,1994)
	keep if default_serial == 1
	drop if return_real==.
	bys bond_ID y: egen countmonths = count(return_real)
	drop if countmonths <= 11
	drop countmonths

	* drop outliers
	drop if country == "Romania" & y == 1952
	drop if country == "Romania" & y == 1959

	tsset bond_ID tm

replace return_real = (1+return_real)*(1-0.41)-1  if currency == "USD" & tm == ym(1934,1)


* Aggregate data 
	collapse (mean) return_real	[aw = amount_issUSD], by(y m tm)
	
* Label asset class	
	gen asset = "External, sovereign bonds"
	gen bond_ID = 0


* Append alternative asset classes 
*S&P 500
	append using "other_asset_classes_returns", keep( returns_SP_real tm y m)
	replace bond_ID = 101 if bond_ID == .
	replace return_real = returns_SP_real if bond_ID == 101
	replace asset = "S&P Return Index" if bond_ID == 101
	drop returns_SP_real
	
* Append alternative asset classes 
*US Bonds
	append using "other_asset_classes_returns", keep(returns_USTreasury_real tm y m)
	replace bond_ID = 105 if bond_ID == .
	replace return_real = returns_USTreasury_real if bond_ID == 105
	replace asset = "Treasury Bonds Return US" if bond_ID == 105
	drop returns_USTreasury_real

* Labeling: Dummy for crises 	
	gen dum_cris = .
	replace dum_cris = 0 if y == 1907 & m == 10
	replace dum_cris = 0 if y == 1929 & m == 10 
	replace dum_cris = 0 if y == 2008 & m == 9 
 
* Labeling: asset names 
	gen c1=""
	replace c1 = "returns_SP_real" 	if asset == "S&P Return Index"
	replace c1 = "EXTSOV" 			if asset == "External, sovereign bonds"
	replace c1 = "returns_USTreasury_real" 		if asset == "Treasury Bonds Return US"

* Labeling: Dummy for crises 	
	gen event = 1 if inrange(tm,ym(1904,10), ym(1910,10))
	replace event = 2 if inrange(tm,ym(1926,10), ym(1937,10))
	replace event = 3 if inrange(tm,ym(2005,9), ym(2011,9))

* Compute cumulative returns 
	global lag = 36

	forvalues t = 1/$lag {
		bys asset (tm): replace dum_cris = `t' if dum_cris[_n-`t'] == 0
	}
	forvalues t = 1/$lag {
		bys asset (tm): replace dum_cris = -`t' if dum_cris[_n+`t'] == 0
	}


	forvalues t = 37/100 {
		bys asset (tm): replace dum_cris = `t' if dum_cris[_n-`t'] == 0 & event==2
	}

	drop if dum_cris == .

	
* Devaluation
	sort bond_ID tm
	by bond_ID: replace return_real = (1+return_real)*(1-0.41)-1  if c1 == "returns_SP_real" & tm == ym(1934,01)
	by bond_ID: replace return_real = (1+return_real)*(1-0.41)-1  if c1 == "returns_USTreasury_real" & tm == ym(1934,01)


*rolling mean/cumret
	gen cumret = 100 if dum_cris == -1
	sort bond_ID event tm dum_cris
	by bond_ID event: replace cumret = cumret[_n-1] * (1+return_real[_n]) if _n > $lag
	gen neg = -dum_cris
	sort bond_ID event neg
	by bond_ID event: replace cumret = cumret[_n-1] * (1/ (1+return_real[_n-1]) ) if _n > ($lag + 2)
	drop neg
	
	

* Produce and export figure 15 panel A: 1907 New York Panic	
* ==========================================================
	sum cumret if event == 1
	gen max = r(max) if dum_cris == -1|dum_cris == 0|dum_cris == -2 & event == 1
	gen min = r(min) if dum_cris == -1|dum_cris == 0|dum_cris == -2 & event == 1
			
	replace min = 50 if dum_cris == -1| dum_cris == 0|dum_cris == -2 & event == 1
	replace max = 150 if dum_cris == -1|dum_cris == 0|dum_cris == -2 & event == 1
	
	
	#delimit;		
	twoway
				(area min  dum_cris if event==1& c1=="returns_SP_real" & min>=50,  cmissing(no) color(gs15)) 
				(area  max dum_cris if event==1& c1=="returns_SP_real"& min>=50, cmissing(no) color(gs15)) 
				(line cumret dum_cris if event==1& c1=="returns_SP_real"& min>=50, lpattern(longdash) lwidth(medthick) color(gs0)) 
				(line  cumret dum_cris if  event==1& c1=="EXTSOV" & min>=50, lwidth(medthick) color(maroon)) 
				(line  cumret dum_cris if  event==1& c1=="returns_USTreasury_real" & min>=50, lpattern(shortdash) lwidth(medthick) color(gs5)) 
				, ylabel(50 75 100 125 150,nogrid angle(0) format(%4.0f) labsize(medlarge)) 
				graphregion(color(white)) 
				ytitle("Cumulative total return index", size(medlarge) margin(right)) 
				xtitle("Years around financial crisis", size(medlarge)) 
				xlabel(-36 "-3" -24 "-2" -12 "-1" 0 "t" 12 "1" 24 "2" 36 "3" , labsize(medlarge)) 
				legend( col(2) size(medlarge) 
				symysize(*.5)  symxsize(*.5) rowgap(*0.01) colgap(*0.01) 
				region(color(white)) 
				label(5 "US Treasuries (10 year)") 
				label(3 "US equities")  
				label(4 "External sovereign bonds")  
				order(3 5 4))


	;
	#delimit cr
	drop min max

	graph export "Figure15_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure15_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure15_PanelA.tif", as(tif) replace
	
* Produce and export figure 15 panel B: 1930 Great Depression
* ==========================================================
	sum cumret 
	gen max=r(max) if dum_cris==-1|dum_cris==0|dum_cris==-2 
	gen min=r(min) if dum_cris==-1|dum_cris==0|dum_cris==-2 
	replace min=50 if dum_cris==-1|dum_cris==0|dum_cris==-2 
	replace max=150 if dum_cris==-1|dum_cris==0|dum_cris==-2 

	#delimit;
	twoway	
				(area min dum_cris if event == 2 & c1 =="returns_SP_real" & dum_cris <=36,  cmissing(no) color(gs15))
				(area max dum_cris if event == 2 & c1 =="returns_SP_real" & dum_cris <=36, cmissing(no) color(gs15)) 
				(line cumret dum_cris if event == 2 & c1 =="returns_SP_real" & dum_cris <=36, lpattern(longdash)  lwidth(medthick) color(gs0)) 
				(line cumret dum_cris if  event == 2 & c1 =="EXTSOV" & dum_cris <=36, lwidth(medthick) color(maroon)) 
				(line cumret dum_cris if  event == 2 & c1 =="returns_USTreasury_real" & dum_cris <=36 , lpattern(shortdash) lwidth(medthick) color(gs5))
				, ylabel(0 25 50 75 100 125 150,nogrid angle(0) format(%4.0f) labsize(medlarge))
				graphregion(color(white)) 
				xsize(6)
				ytitle("Cumulative total return index", size(medlarge) margin(right)) 
				xtitle("Years around financial crisis", size(medlarge)) 
				xlabel(-36 "-3" -24 "-2" -12 "-1" 0 "t" 12 "1" 24 "2" 36 "3" , labsize(medlarge)) 
				legend( col(2) size(medlarge) 
				symysize(*.5)  symxsize(*.5) rowgap(*0.01) colgap(*0.01) 
				region(color(white)) 
				label(5 "US Treasuries (10 year)") 
				label(3 "US equities")  
				label(4 "External sovereign bonds")  
				order(3 5 4))
		
			
	;
	#delimit cr
			
	graph export "Figure15_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure15_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure15_PanelB.tif", as(tif) replace
	drop min max		

* Produce and export figure 15 panel C: Great Recession 2007/2008
* ==============================================================================
	gen max=150 if (dum_cris==-1|dum_cris==0|dum_cris==-2) 
	gen min=50 if (dum_cris==-1|dum_cris==0|dum_cris==-2 ) 

	#delimit;	
	twoway
				(area min  dum_cris if event==3& c1=="returns_SP_real" ,  cmissing(no) color(gs15)) 
				(area  max dum_cris if event==3& c1=="returns_SP_real", cmissing(no) color(gs15)) 
				(line cumret  dum_cris if event==3& c1=="returns_SP_real", lpattern(longdash)  lwidth(medthick) color(gs0)) 
				(line  cumret dum_cris if  event==3& c1=="EXTSOV" , lwidth(medthick) color(maroon)) 
				(line  cumret dum_cris if  event==3& c1=="returns_USTreasury_real" , lpattern(shortdash) lwidth(medthick) color(gs5)) 
				, ylabel(50 75 100 125 150,nogrid angle(0) format(%4.0f) labsize(medlarge)) 
				graphregion(color(white)) 
				ytitle("Cumulative total return index", size(medlarge) margin(right)) 
				xtitle("Years around financial crisis", size(medlarge)) 
				xlabel(-36 "-3" -24 "-2" -12 "-1" 0 "t" 12 "1" 24 "2" 36 "3" , labsize(medlarge)) 
				legend( col(2) size(medlarge) 
				symysize(*.5)  symxsize(*.5) rowgap(*0.01) colgap(*0.01) 
				region(color(white)) 
				label(5 "US Treasuries (10 year)") 
				label(3 "US equities")  
				label(4 "External sovereign bonds")  
				order(3 5 4))
	;
	#delimit cr

	drop min max
	graph export "Figure15_PanelC.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure15_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure15_PanelB.tif", as(tif) replace


	
*******************************************************************
****************** Part III: Tables *************************************
*******************************************************************

* =================================================================
* Table 2
* =================================================================
* Load main dataset
	use returns_prices_1820_2016.dta, clear

* Restrict to active bonds
	drop if price==.

* Compute maturity
	gen maturityp = maturityy - issuey +1 if maturityy!=. & issuey!=.
	
* Compute dummy for eras
	gen decade=1 if inrange(y,1815,1869)
	replace decade=2 if inrange(y,1870,1913)
	replace decade=3 if inrange(y,1914,1945)
	replace decade=4 if inrange(y,1946,1980)
	replace decade=5 if inrange(y,1989,2016)

	gen bond_IDS=string(bond_ID)

* Calculate bond pricing and countries by era, loop through eras. 

foreach i in 1 2 3 4 5 6 {								 // loop through decades
local r=1

 if `i' == 6 {
   replace decade=6 									// decade = 6 --> entire period
  }
   
   
mat M`i' = J(9,1,.)

* Number of countries active in market
	quietly tab country if decade==`i'
	mat M`i'[`r',1] = r(r) 
	local r = `r'+1

* Share of country active in bond market to total independent countries
	preserve
	keep if decade == `i'
	quietly tab country
	local nc = r(r)
	bys y: keep if _n == 1
	quietly sum NRcIndependent 
	local m = r(max) 
	mat M`i'[`r',1] = real(string(( `nc'/`m')*100,"%4.1f"))
	restore
	local r = `r'+1

* Number price per decade
	quietly count if price ~=. & decade ==`i'
	mat M`i'[`r',1] = r(N) 
	local r = `r'+1

* Number of active bonds 
	quietly tab bond_ID if decade == `i'
	mat M`i'[`r',1] =  r(r)
	local r = `r'+1

* Number of active GBP bonds 
	quietly tab bond_ID if decade == `i' & currency =="GBP"
	mat M`i'[`r',1] = r(r)
	local r = `r'+1

* Number of active USD bonds 
	quietly tab bond_ID if decade == `i' &  currency =="USD"
	mat M`i'[`r',1] = r(r)
	local r = `r'+1

 *Average Maturity
	preserve
	keep if decade == `i'
	bys bond_ID: keep if _n == 1
	quietly sum maturityp 
	mat M`i'[`r',1] = real(string(r(mean),"%4.2f")) 
	restore
	local r = `r'+1

*Average Coupon of active bond within each decade
	preserve
	keep if decade == `i'
	bys bond_ID: keep if _n == 1 
	quietly sum coupon 
	mat M`i'[`r',1] = real(string(r(mean),"%4.2f"))
	restore
	local r = `r'+1

*Average Amount issued of active bond within each decade
	preserve
	keep if decade == `i'
	bys bond_ID: keep if _n == 1 
	quietly sum amount_issUSD 
	mat M`i'[`r',1] = real(string(r(mean),"%4.2f")) 
	restore

}

* Put excel 	
* ---------------------------------------------------------
		
* Define row names 	
	#delimit;
	local row `" "Number of countries covered" 
				"Share of countries covered"
				"Pricing observations"
				"Number of active bonds" 
				"… issued in British pounds"
				"… issued in US dollars"
				"Average maturity of bonds"
				"Average coupon (nominal)"
				"Average amount issued (nominal, in m USD)" "'
	;
	#delimit cr
	
	
* Open output file  	
	putexcel set "tables.xlsx", sheet("Table 2")  modify

* Export row names 	
	local i = 4
	foreach d of local row{ 								 
	putexcel B`i' = ("`d'")
	local i = `i' +1
	}
	
* Export stats 														
	putexcel (C4)= matrix(M6)					
	putexcel (E4)= matrix(M1)
	putexcel (F4)= matrix(M2)
	putexcel (G4)= matrix(M3)
	putexcel (H4)= matrix(M4)
	putexcel (I4)= matrix(M5)
	
	
* Format excel cells	
	$format1	
	putexcel C3 = ("Total sample")
	putexcel F2 = ("By era")
	putexcel E3 = ("1815-1869"), txtwrap $format2 hcenter
	putexcel F3 = ("1870-1913"), txtwrap $format2 hcenter
	putexcel G3 = ("1914-1945"), txtwrap $format2 hcenter
	putexcel H3 = ("1946-1980"), txtwrap $format2 hcenter
	putexcel I3 = ("1989-2016"), txtwrap $format2 hcenter


	putexcel B2:H2, border(top,double) $format2 overwritefmt 
	putexcel B4:H4, border(top) $format2 overwritefmt
	putexcel B13:H13, border(top,double) $format2 overwritefmt

	$format1	
	mata													
	b = xl()
	b.load_book("tables.xlsx")
	b.set_sheet("Table 2")
	b.set_column_width(2,2, 35)
	b.set_column_width(3,3, 12)
	b.set_column_width(4,8, 10)
	b.close_book()
	end
	
	putexcel C3, bold hcenter
	putexcel F2, bold hcenter


* =================================================================
* Table 3 : Returns on a global portfolio of external sovereign bonds, 1815-2016
* =================================================================

* Load time series returns data 
	use loadDataTS, clear
	
* Compute monthly summary stats 
 
	gen t=.
	foreach var of varlist return_real return_nom { 				// compute geometric returns
	gen lret = log(`var'+1)
	egen sumlret = total(lret)
	gen exp`var' = exp(sumlret)
	egen county`var' = count(lret)
	drop sumlret lret
	}

	gen geomR = ((expreturn_real)^(1/countyreturn_real)-1)
	gen geomN = ((expreturn_nom)^(1/countyreturn_nom)-1)
	drop countyreturn_nom countyreturn_real expreturn_real expreturn_nom

	egen SDr = sd(return_real) 							// compute standareturn_nom deviation 
	egen SDn = sd(return_nom)
	
	egen meanExcess = mean(excessBonds)
	gen sharpeMonthly = meanExcess/SDr

	mat M0 = J(1,10,.)							// put monthly summary stats in matrix
	local c = 1
	foreach var of varlist return_real geomR SDr t return_nom geomN SDn t excessBonds sharpeMonthly{
	quietly sum `var' 
	mat M0[1,`c'] = real(string(r(mean)*100,"%4.2f"))  
	local c =`c'+1
	}


* Compute yearly summary stats 
collapse ryreal ry ryrealExcess, by(y)

* Create dummy for eras
	gen decade=.
	replace decade=1 if inrange(y,1815,1869)
	replace decade=2 if inrange(y,1870,1914)
	replace decade=3 if inrange(y,1915,1945)
	replace decade=4 if inrange(y,1946,1973)
	replace decade=5 if inrange(y,1995,2016)
	
* Create dummy for World War I and World War II 
			
	gen WW=.
	replace WW=1 if inrange(y,1914,1918)	
	replace WW=1 if inrange(y,1938,1945)	

* Compute yearly summary stats 
	
	gen t=.
				
	foreach i in 1 2 3 4 5 6 7{							 // loop through decades
	local r=1
	preserve 

	 if `i' == 6 {
	   replace decade=6 								// decade = 6 --> entire period
	  }
	   

	 if `i' == 7 {										// decade = 7 --> entire period without World Wars
	   replace decade=7
	   drop if WW==1
	  }
	  
	keep if decade==`i'


	foreach var of varlist ry ryreal{				// compute geometric returns
	gen lret=log(`var'+1)
	egen sumlret=total(lret)
	gen exp`var'=exp(sumlret)
	egen county`var'=count(lret)
	drop sumlret lret
	}

	gen geomR=((expryreal)^(1/countyryreal)-1)
	gen geomN=((expry)^(1/countyry)-1)
	drop expry countyry expryreal countyryreal

	egen SDr=sd(ryreal)								// compute standareturn_nom deviation 
	egen SDn=sd(ry)

	egen sd_sharpeSafeSynthetic = sd(ryrealExcess)		// compute sharpe ratio
	egen excessSafeSynthetic_average = mean(ryrealExcess)
	gen sharpeSafeSynthetic = excessSafeSynthetic_average / sd_sharpeSafeSynthetic
	drop sd_sharpeSafeSynthetic
		

	mat M`i' = J(1,10,.)								// put yearly summary stats in matrix
	local c = 1
	foreach var of varlist ryreal geomR SDr t ry geomN SDn t excessSafeSynthetic_average{
	quietly sum `var' 
	mat M`i'[1,`c'] = real(string(r(mean)*100,"%4.2f"))  
	local c =`c'+1
	}

	quietly sum sharpeSafeSynthetic 
	mat M`i'[1,10] = real(string(r(mean),"%4.2f"))  


	restore
	}

* Put excel 	
* ---------------------------------------------------------
	
* Define row names 	
	#delimit;
	local row `" "Full sample, 1815-2016, yearly" 
				"	… without world wars"
				" "
				"… without world wars" 
				"	1815-1869"
				"	1870-1914"
				"	1915-1945"
				"	1946-1973"
				"	1995-2016"
				" "
				"Monthly returns, full sample" "'
	;
	#delimit cr
	
* Open output file 	
	putexcel set "tables.xlsx", sheet("Table 3")  modify    // open output file
	
	
* Export row names 
	local i = 6
	foreach d of local row{ 							
	putexcel B`i'=("`d'")
	local i = `i' +1
	}
		
* Export stats 
	putexcel (C10) = matrix(M1) 						
	putexcel (C11) = matrix(M2) 		
	putexcel (C12) = matrix(M3) 	
	putexcel (C13) = matrix(M4) 	
	putexcel (C14) = matrix(M5)		
	putexcel (C6)  = matrix(M6)		
	putexcel (C7)  = matrix(M7)			
	putexcel (C16) = matrix(M0)		
	
* Format excel cells
	$format1
	putexcel C4:E4 = "Real Return", merge hcenter border(bottom) $format2
	putexcel C5 = "Arithm. Mean",  txtwrap hcenter top $format2
	putexcel D5 = "Geom. Mean",  txtwrap hcenter top  $format2
	putexcel E5 = "SD"
	putexcel G5 = "Arithm. Mean"
	putexcel H5 = "Geom. Mean"
	putexcel I5 = "SD"
	putexcel G4:I4 = "Nominal Return", merge hcenter border(bottom) $format2
	putexcel K4:K5 = "Above US/UK gov. Bonds", merge txtwrap hcenter top $format2
	putexcel L4:L5 = "Sharpe ratio", merge txtwrap hcenter top $format2
	putexcel B3:L3, border(bottom,double)$format2 overwritefmt
	putexcel B5:L5, txtwrap hcenter top border(bottom) $format2 overwritefmt
	putexcel B16:L16, border(bottom,double) $format2 overwritefmt
	
	
	mata													
	b = xl()
	b.load_book("tables.xlsx")
	b.set_sheet("Table 3")
	b.set_column_width(2,2,30) 		//make  column wbond_IDest
	b.set_column_width(6,6,1) 		//make  column less wbond_IDe
	b.set_column_width(10,10,1) 	//make  column less wbond_IDe
	b.set_row_height(5,5,50) 	//make  row hight larger
	b.close_book()
	end


* =================================================================
* Table 4
* =================================================================

* Load time series returns data 
	use loadDataTS, clear

* -------------------------------------------------------------------------
	gen bond_ID = 1
	tsset bond_ID tm 

	forvalues i=1/4{
	sort bond_ID tm 
	local ii = ((3*`i')-1)
	gen double lret1 = log(return_real+1)
	rangestat (sum) lret1 (count) lret1 (sd)lret1, interval(tm -`ii' 0) by(bond_ID)
	gen expreturn_real`i' = exp(lret1_sum)
	gen arithQ_`i' = exp(lret1_sum) - 1 if lret1_count == `ii'+1
	drop lret1 lret1_count lret1_sd exp* lret1_sum
	}


	gen arithQ = arithQ_1
	replace arithQ_1 =. if m ~= 3
	replace arithQ_2 =. if m ~= 6
	replace arithQ_3 =. if m ~= 9
	replace arithQ_4 =. if m ~= 12

* Compute quarterly geometric returns 

	gen geom2 = (1+arithQ_2)^(1/((2*3)/3))-1
	gen geom3 = (1+arithQ_3)^(1/((3*3)/3))-1
	gen geom4 = (1+arithQ_4)^(1/((4*3)/3))-1

* Compute quarterly average returns 
	gen arithm2 = arithQ_2/2
	gen arithm3 = arithQ_3/3
	gen arithm4 = arithQ_4/4

* Compute summary stats and put everything in a matrix
	local i = 1
	mat M`i' = J(6,4,.)
	local stats "mean p75 p50 p25"
	matrix rownames M`i' = `stats'
	local c = 1 


	foreach var of varlist arithQ_1 arithQ_2 arithQ_3 arithQ_4 {
	local r = 1
	foreach stat of local stats  {
	qui sum `var' if `var' ~=. , detail
	mat M`i'[`r',`c'] = real(string(r(`stat')*100,"%4.1f")) 
	local r = `r'+1
	}
	local c = `c'+1
	}

	local c = 1
	foreach var of varlist arithQ_1 arithm2 arithm3 arithm4{
	qui sum `var' if `var' ~=. , detail
	mat M`i'[5,`c'] = real(string(r(mean)*100,"%4.1f")) 
	local c = `c'+1
	}


	local c = 1
	foreach var of varlist arithQ_1 geom2 geom3 geom4{
	qui sum `var' if `var' ~=. , detail
	mat M`i'[6,`c'] = real(string(r(mean)*100,"%4.1f"))
	local c = `c'+1
	}

	

* Compute holding period summary stats 
* -------------------------------------------------------------------------

* Compute yearly returns  
	collapse ryreal ry, by(y)	

* Compute annual holding period returns 
	gen bond_ID = 1
	tsset bond_ID y 

	foreach var in ryreal{

	forvalues i = 1/10{
	sort bond_ID  y 
	local ii = ((1*`i')-1)
	gen double lret1 = log(`var'+1)
	rangestat (sum) lret1 (count) lret1 (sd)lret1, interval(y -`ii' 0) by(bond_ID)
	gen exp`var'`i' 	= exp(lret1_sum)
	gen arith_`var'`i' 	= exp(lret1_sum) - 1 				if lret1_count == `ii'+1
	gen geom_`var'`i' 	= exp`var'`i'^(1/(lret1_count)) - 1 if lret1_count == `ii'+1
	gen  sd_`var'`i'	= lret1_sd 							if lret1_count == `ii'+1
	drop lret1 lret1_sum lret1_count lret1_sd exp`var'`i'
	replace geom_`var'`i'	= geom_`var'`i'
	replace arith_`var'`i'	= arith_`var'`i'
	gen sdAnnual_`var'`i'	= (sd_`var'`i'*sqrt(12))
	}
	}



* Compute annual arithmetic and geometric returns for different holding period returns 

	foreach var in arith_ryreal1{

	forvalues i = 1/10{
	sort bond_ID y 
	local ii = ((1*`i')-1)
	gen double lret1 = log(`var'+1)
	rangestat (sum) lret1 (mean) arith_ryreal1 (median) arith_ryreal1 (sd) arith_ryreal1 (count) arith_ryreal1, interval(y -`ii' 0) by(bond_ID)
	gen return_nommean`i' = arith_ryreal1_mean 						   if arith_ryreal1_count ==`ii'+1
	gen return_nommedian`i' = arith_ryreal1_median  				   if arith_ryreal1_count ==`ii'+1
	gen return_nommeanSD`i' = arith_ryreal1_sd 						   if arith_ryreal1_count ==`ii'+1
	gen return_nommeanNR`i' = arith_ryreal1_count   				   if arith_ryreal1_count ==`ii'+1
	gen exp`var'`i' = exp(lret1_sum)  						   if arith_ryreal1_count ==`ii'+1
	gen geom`i' = (exp`var'`i'^(1/(arith_ryreal1_count)) - 1)  if arith_ryreal1_count ==`ii'+1
	drop arith_ryreal1_sd arith_ryreal1_median arith_ryreal1_mean arith_ryreal1_count exp`var'`i' lret1 lret1_sum
	}

	}


* Compute summary stats and put everything in a matrix 
	mat M2 = J(6,10,.)
	matrix rownames M2 = `stats'
	local c = 1 

	foreach var of varlist arith_ryreal* {
	local r = 1
	foreach stat of local stats  {
	qui sum `var' if `var' ~=. , detail
	mat M2[`r',`c'] = real(string(r(`stat')*100,"%4.1f")) 
	local r = `r'+1
	}
	local c = `c'+1
	}
	local c = 1 
	foreach var of varlist return_nommean1 return_nommean2 return_nommean3 ///
	return_nommean4 return_nommean5 return_nommean6 return_nommean7 ///
	return_nommean8 return_nommean9 return_nommean10 {
	qui sum `var'
	mat M2[5,`c'] = real(string(r(mean) *100,"%4.1f")) 
	local c = `c'+1
	}

	local c = 1 
	foreach var of varlist geom_ryreal* {
	qui sum `var'
	mat M2[6,`c'] = real(string(r(mean) *100,"%4.1f"))
	local c = `c'+1
	} 
		
	
* Put excel 
* ----------------------------------------------------------------
* Define row names 	
	#delimit;
	local row `" "Mean (arithmetic,cumulative)"
				"   p75"
				"   Median" 
				"	p25"
				"Mean (arithmetic, avg. per period)"
				"Mean (geometric,annualized)" "'
	;
	#delimit cr
	
* Open output file	
	putexcel set "tables.xlsx", sheet("Table 4")  modify 
	
* Export row names 
	local i = 8
	foreach d of local row{ 							
	putexcel B`i'=("`d'") , txtwrap
	local i = `i' +1
	}
		
	
* Export stats 
	putexcel (C8) = matrix(M1)
	putexcel (H8) = matrix(M2) 

* Format excel cells
	$format1
	putexcel C6:F6 = "Quarterly (Q1-Q4)", merge hcenter $format2
	putexcel H6:Q6 = "Yearly (years 1-10)", merge left  $format2
	
	local i = 1
	foreach col in C D E F {
	putexcel `col'7 = "Q`i'", hcenter $format2
	local i = `i' +1
	}
	
	local i = 1
	foreach col in H I J K L M N O P Q {
	putexcel `col'7 = "`i'",  hcenter $format2
	local i = `i' +1
	}
	
	putexcel B5:Q5, border(bottom,double) $format2 overwritefmt
	putexcel B14:Q14, border(top,double) $format2 overwritefmt	
	putexcel B7:Q7, txtwrap hcenter top border(bottom) $format2 overwritefmt
	
	
	mata													
	b = xl()
	b.load_book("tables.xlsx")
	b.set_sheet("Table 4")
	b.set_column_width(7,7,1) 		//make  column less wbond_IDe
	b.set_column_width(2,2,30) 	//make  row hight larger
	b.close_book()
	end

	
	
* =================================================================
* Table 5: Bond returns by country group and era
* =================================================================

* Load returns data 
	use loadData, clear
	
* Compute return portfolios by default history: serial defaulters vs. non-serial defaulters	
	#delimit ;	
	collapse (mean) return_nom return_real excessbillsredeal excessreturn_real 
	[aw = amount_issUSD], by(default_serial y m tm)
	;
	#delimit cr	
	
* Restrict 	
	drop if return_real==.
	bys default_serial y: egen countmonths = count(m)
	drop if countmonths != 12
	drop countmonths

* Compute annual returns 
	foreach var in 	return_nom  return_real excessbillsredeal excessreturn_real{
	bys default_serial y: egen lryreal = total( ln( `var'+1 ) ), missing
	gen `var'y = exp( lryreal ) - 1 
	drop lryreal	
	}	
		

collapse (mean) return_nomy return_realy excessbillsredealy excessreturn_realy, by(default_serial y)

	
* Add world risk-free rate 	
	append using "other_asset_classes_returns_yearly.dta"	, keep(rySafe_real y) gen(marker)
	replace default_serial = 2 if marker == 1
	replace return_realy = rySafe_real if marker == 1
	drop marker rySafe_real
		
* Compute geometric return 		
	foreach var of varlist return_realy{
	bys default_serial: gen lret = log(`var'+1)
	bys default_serial: egen sumlret = total(lret)
	gen exp`var' = exp(sumlret)
	bys default_serial: egen county`var' = count(lret)
	drop sumlret lret
	}

	gen geomR=((expreturn_realy)^(1/countyreturn_realy)-1)
	bys default_serial: egen SDr = sd(return_realy)
	drop expreturn_realy countyreturn_realy


* Compute dummy for decades
	gen decade =.
	replace decade=1 if inrange(y,1815,1869)
	replace decade=2 if inrange(y,1870,1913)
	replace decade=3 if inrange(y,1920,1937)
	replace decade=4 if inrange(y,1995,2016)
	
* Compute geometric returns by decades

	foreach var of varlist return_realy {
	bys default_serial decade: gen lret = log(`var'+1)
	bys default_serial decade: egen sumlret = total(lret)
	gen exp=exp(sumlret)
	bys default_serial decade: egen county = count(lret)
	bys default_serial decade: egen SDdecade`var' = sd(`var')
	gen geomreturn_nomecade`var' = ((exp)^(1/county)-1)
	capture drop lret sumlret exp county
	}

* Compute summary stats and put everything in a matrix

	foreach e in 1 2 3 4 5{
	if `e' == 5 replace decade =`e'
	mat M`e' = J(3,5,.)
	local r = 1 
	local stats "N mean p50 mean mean"

	foreach i in 1 0 2{
	local c = 1 
	quietly sum return_realy if decade ==`e' & default_serial ==`i'
	mat M`e'[`r',`c'] = (real(string(100*r(mean),"%4.1f")))
	local c = `c'+1
	quietly sum geomreturn_nomecadereturn_realy if decade ==`e' & default_serial ==`i'
	mat M`e'[`r',`c'] = (real(string(100*r(mean),"%4.1f"))) 
	local c = `c'+1
	quietly sum SDdecadereturn_realy if decade ==`e' & default_serial ==`i'
	mat M`e'[`r',`c'] = (real(string(100*r(mean),"%4.1f")))
	local c = `c'+1
	quietly sum excessreturn_realy if decade ==`e' & default_serial ==`i'
	mat M`e'[`r',`c'] = (real(string(100*r(mean),"%4.1f")))

	if `e'==5 {
		sum return_realy if default_serial ==`i'
		mat M`e'[`r',1] = (real(string(100*r(mean),"%4.1f"))) 
		sum excessreturn_realy if default_serial ==`i'
		mat M`e'[`r',4] = (real(string(100*r(mean),"%4.1f")))
		preserve
			keep if decade == `e'
			bys default_serial: keep if _n == 1
			sum geomR if default_serial == `i'
			mat M`e'[`r',2] = (real(string(100*r(mean),"%4.1f")))
			sum SDr if  default_serial ==`i'
			mat M`e'[`r',3] = (real(string(100*r(mean),"%4.1f"))) 
		restore
		}
	local r = `r'+1	
	}
	}	

	
* Put excel 	
* ---------------------------------------------------------
	
* Define row names 	
	#delimit;
	local row `" "Total Sample (1815-2016)" 
				" Serial defaulters"
				" Other countries with ext. bonds"
				" UK/US government bonds" 
				" "
				"Early 19th Century (1815-1869)"
				" Serial defaulters"
				" Other countries with ext. bonds"
				" UK government bonds" 
				" "
				"Pre-WW1 (1870-1913)"
				" Serial defaulters"
				" Other countries with ext. bonds"
				" UK government bonds" 
				" "
				"Interwar (1920-1937)"
				" Serial defaulters"
				" Other countries with ext. bonds"
				" US government bonds"
				" "
				"Today (1995-2016)"
				" Serial defaulters"
				" Other countries with ext. bonds"
				" US government bonds" "'
		
	;
	#delimit cr
	
* Open output file 
	putexcel set "tables.xlsx", sheet("Table 5")  modify   
	
* Export row names 
	local i = 6
	foreach d of local row{ 							
	putexcel B`i' = ("`d'")
	local i = `i' +1
	}
	
* Export stats 
	putexcel (C12) = matrix(M1) 
	putexcel (C17) = matrix(M2) 
	putexcel (C22) = matrix(M3) 
	putexcel (C27) = matrix(M4) 
	putexcel (C7)  = matrix(M5) 
	
* Format excel cells
	$format1
	putexcel C5 = "Arithmetic mean (annual)"
	putexcel D5 = "Geometric mean (annual)"
	putexcel E5 = "SD"
	putexcel F5 = "Excess return (mean, above UK/US bonds)"
	
	putexcel B4:F4, border(bottom,double) $format2 overwritefmt
	putexcel B30:F30, border(top,double) $format2 overwritefmt	
	putexcel B5:F5, txtwrap hcenter top border(bottom) $format2 overwritefmt
	
	mata													
	b = xl()
	b.load_book("tables.xlsx")
	b.set_sheet("Table 5")
	b.set_column_width(2,2,30) //make row column wbond_IDest
	b.set_column_width(3,4,10) //make row column wbond_IDest
	b.set_column_width(6,6,13) //make row column wbond_IDest
	b.set_row_height(5,5,50) 	//make  row hight larger
	b.close_book()
	end

	
* =================================================================
* Table 6: Sovereign bond returns by country
* =================================================================

* Load country return data 
	use loadDatac, clear
		
* Compute standareturn_nom deviation of annual arithmetic returns
	bys country (y): egen sd = sd(return_realy)

* Compute sharpe ratio over bills synthetic
	by country: egen sd_sharpebillsSynthetic = sd(excessbillsredealy)
	by country: egen excessbillsSynthetic_average = mean(excessbillsredealy)
	gen sharpebillsSynthetic = excessbillsSynthetic_average / sd_sharpebillsSynthetic


* Compute sharpe ratio over bonds synthetic
	by country: egen sd_sharpebondsSynthetic = sd(excessbillsredealy)
	by country: egen excessbondsSynthetic_average = mean(excessbillsredealy)
	gen sharpebondsSynthetic = excessbillsSynthetic_average / sd_sharpebillsSynthetic

* Compute geometric returns entire sample
	foreach var of varlist return_realy{
	gen lret = log(`var'+1)
	by country: egen sumlret = total(lret)
	gen exp`var' = exp(sumlret)
	bys country: egen county`var' = count(lret)
	drop sumlret lret
	}

	gen geomR = ((expreturn_realy)^(1/countyreturn_realy)-1)
	drop expreturn_realy countyreturn_realy

* Compute geometric returns by subsample
	gen sample = 1 if y<1980

	foreach var of varlist return_realy{
	gen lret = log(`var'+1)
	bys country sample: egen sumlret = total(lret)
	gen exp`var' = exp(sumlret)
	bys country sample: egen county`var' = count(lret)
	drop sumlret lret
	}
	
* Add historical geometric return
	preserve
	keep if sample == 1
	gen geomRhist = ((expreturn_realy^(1/countyreturn_realy))-1)
	bys country: keep if _n == 1
	keep country geomRhist
	tempfile temp
	save "`temp'", replace
	restore

	merge m:1 country using `temp'
	drop _merge 
	
* Add modern geometric return
	preserve
	keep if sample == .
	gen geomRmodern=((expreturn_realy)^(1/countyreturn_realy)-1)
	bys country: keep if _n == 1
	keep country geomRmodern
	tempfile temp
	save "`temp'", replace
	restore

	merge m:1 country using `temp'
	drop _merge 
	
* Add modern + historical arithmetic return	
	bys country sample: egen Meanhist = mean(return_realy)
	replace Meanhist = . if sample == .

	bys country sample: egen MeanModern = mean(return_realy)
	replace MeanModern = . if sample == 1
	
* Add skewness
	bys country: egen skewness = skew(return_realy)

* Drop outlier 
	drop if country == "Romania" & y == 1952
	drop if country == "Romania" & y == 1959

	
* Collapse 	
	bys country y: keep if _n == 1
	

	#delimit ;	
	collapse (mean) 
	sharpebillsSynthetic
	sharpebondsSynthetic
	excessbillsSynthetic_average
	excessbondsSynthetic_average
	return_realy
	sd
	sd_sharpebillsSynthetic
	sd_sharpebondsSynthetic
	geomR geomRhist geomRmodern
	MeanModern Meanhist
	skewness
	(count) count = return_realy
	, by(country)
	;
	#delimit cr

* Restrict 
	drop if count < 10
	gsort -return_realy
	gen n = _n
	order country return_realy MeanModern Meanhist

* Only returns data for modern period
	levelsof country if Meanhist == . & MeanModern !=. , local(c1) 
	foreach c of local c1 {
	replace country = country + "°" if country == "`c'"
	}
	
* Only returns data for historical period	
	levelsof country if Meanhist !=. & MeanModern ==. , local(c1) 
	foreach c of local c1 {
	replace country = country + "*" if country == "`c'"
	}
	
* Add info to country names
	replace country = "France (interwar)*"   if country == "France*"
	replace country = "Sweden (pre-WW2)*"    if country == "Sweden*"
	replace country = "Finland (interwar)*"  if country == "Finland*"
	replace country = "Serbia/Yugoslavia"    if country == "Yugoslavia"
	
* Put everything in a matrix	
	qui count
	mat M = J(`r(N)',14,.)

	forvalues i = 1/`r(N)'{
	mat M[`i',1]  = (real(string(return_realy[`i']*100,"%4.1f")))
	mat M[`i',2]  = (real(string(Meanhist[`i']*100,"%4.1f")))
	mat M[`i',3]  = (real(string(MeanModern[`i']*100,"%4.1f")))
	mat M[`i',5]  = (real(string(geomR[`i']*100,"%4.1f")))
	mat M[`i',6]  = (real(string(geomRhist[`i']*100,"%4.1f")))
	mat M[`i',7]  = (real(string(geomRmodern[`i']*100,"%4.1f")))
	mat M[`i',9]  = (real(string(excessbondsSynthetic_average[`i']*100,"%4.1f")))
	mat M[`i',10]  = (real(string(sd[`i']*100,"%4.1f")))
	mat M[`i',11]  = (real(string(sharpebondsSynthetic[`i'],"%4.2f")))
	mat M[`i',12]  = (real(string(skewness[`i'],"%4.2f")))
	mat M[`i',13]  = (real(string(count[`i'],"%4.0f")))
	}
		
* Cross-country average		
	mat MA = J(1,14,.)
	gen t = . 
	
	local i = 1
	foreach var of varlist return_realy Meanhist MeanModern t geomR geomRhist ///
	geomRmodern t excessbondsSynthetic_average sd{
	qui sum `var', det
	mat MA[1,`i']  = (real(string(r(mean)*100,"%4.1f")))
	local i = `i' + 1
	}

	local i = 11
	foreach var of varlist sharpebondsSynthetic skewness count{
	qui sum `var', det
	mat MA[1,`i']  = (real(string(r(mean),"%4.2f")))
	local i = `i' + 1
	}
	
* Cross-country median		
	mat MM = J(1,14,.)
	
	local i = 1
	foreach var of varlist return_realy Meanhist MeanModern t geomR geomRhist ///
	geomRmodern t excessbondsSynthetic_average sd{
	qui sum `var', det
	mat MM[1,`i']  = (real(string(r(p50)*100,"%4.1f")))
	local i = `i' + 1
	}	
	
	
	local i = 11
	foreach var of varlist sharpebondsSynthetic skewness count{
	qui sum `var', det
	mat MM[1,`i']  = (real(string(r(p50),"%4.2f")))
	local i = `i' + 1
	}
	
	drop t
	
* Put excel 
* --------------------------------------------------------------------
* Open output file 	
	putexcel set "tables.xlsx", sheet("Table 6")  modify 
	
* Export row names 
	levelsof n, local(n1)
	local i = 4
	foreach n of local n1{
	putexcel B`i' = country[`n']
	local i = `i' + 1
	}
	
* Export stats 	
	putexcel (C4) = matrix(M)
	putexcel (C70) = matrix(MA)
	putexcel (C71) = matrix(MM)
	
	
* Format excel cells
	$format1
	putexcel O2:O3 = "Years with bond returns", merge top txtwrap
	putexcel B2:B3 = "Country", merge top txtwrap
	putexcel C2:E2 = "Real return", merge $format2 hcenter top  border(bottom)
	putexcel G2:I2 = "Geometric return", merge $format2 hcenter top  border(bottom)
	
	foreach col in C G K L M N {
	putexcel `col'3 = "Full sample", $format2  
	}
	
	foreach col in D H {
	putexcel `col'3 = "1820-1973", 
	}
	
	foreach col in E I{
	putexcel `col'3 = "1994-2016", 
	}
	
	putexcel K2 = "Excess return", $format2 hcenter top txtwrap border(bottom)
	putexcel L2 = "SD of real returns", $format2 hcenter top txtwrap border(bottom)
	putexcel M2 = "Sharpe ratio", $format2 hcenter top txtwrap border(bottom)
	putexcel N2 = "Skew. of real return", $format2 hcenter top txtwrap border(bottom)
	putexcel B70 = "Cross-country average", $format2 hcenter top txtwrap
	putexcel B71 = "Cross-country median", $format2 hcenter top txtwrap 
	

	putexcel B1:O1, border(bottom,double) $format2 overwritefmt
	putexcel B72:O72, border(top,double) $format2 overwritefmt	
	putexcel B3:O3, txtwrap hcenter top border(bottom) $format2 overwritefmt
	
	mata													
	b = xl()
	b.load_book("tables.xlsx")
	b.set_sheet("Table 6")
	b.set_column_width(2,2,20) //make row column wbond_IDest
	b.set_column_width(6,6,1) //make row column wbond_IDest
	b.set_column_width(10,10,1) //make row column wbond_IDest
	b.set_row_height(2,2,50) //make row column wbond_IDest
	b.close_book()
	end


	
	

* =====================================================================
* Figure 14 + Table 7: Asset classes across 200 years: Risk and return
* =====================================================================
* Load time series level portfolio returns data
	use loadDataTS, clear

* Prepare data. 
* Again, we will not drop the unmatched data entries given that
* we have no observation for bond returns during 1815-1922,1914, 1974-1994. Yet,
* the alternative asset classes have entries during these years in addition to
* the coverage years for bond returns dataset. Therefore, we keep the dataset
* with the missing entries for bond returns

	rename return_real returns_ExtBonds_real
	rename return_nom returns_ExtBonds_nom
	rename ryrealExcess ryrealExcessExtBonds
	rename excessBonds excessBondsExtBonds
	rename excessBills excessBillsExtBonds
	
	keep y m tm *ExtBonds* 
	
	merge 1:1 y tm m using "other_asset_classes_returns.dta"
	drop _merge
	drop if y<1800
	
	foreach var of varlist returns_*_real{
	local new : subinstr local var "_real" ""
	rename `var' real_`new'
	}
	
	foreach var of varlist returns_*_nom{
	local new : subinstr local var "_nom" ""
	rename `var' nom_`new'
	}
	
	reshape long real_returns_ nom_returns_  excessBonds excessBills ///
		ry ryreal ryrealExcess, i(y tm m) string
	rename _j ticker 
	
* Compute yearly returns 	
	
	bys ticker y: egen lryreal = total( ln( real_returns_+1 ) ), missing
	replace ryreal = exp( lryreal ) - 1 if lryreal ~=.
	drop lryreal

	bys ticker y: egen lryreal = total( ln( nom_returns_+1 ) ), missing
	replace ry = exp( lryreal ) - 1 if  lryreal ~=.
	drop lryreal
	
	bys ticker y: egen lryreal = total( ln( excessBonds+1 ) ), missing
	gen ryexcessBonds = exp( lryreal ) - 1 if lryreal ~=.
	drop lryreal

	bys ticker y: egen lryreal = total( ln( excessBills+1 ) ), missing
	gen ryexcessBills = exp( lryreal ) - 1 if  lryreal ~=.
	drop lryreal
	
	drop excessBills excessBonds
	rename ryexcessBills excessBills
	rename ryexcessBonds excessBonds
	

collapse ryreal ry excess* , by(y ticker)


* Create world-safe asset 	
	
	foreach asset in UKlongbond USTreasury{
	preserve
	keep if  ticker == "`asset'"
	keep if (y<1919  &  ticker == "UKlongbond") | (y>=1919 &  ticker == "USTreasury")
	replace ticker = "safe"
	tempfile temp
	save "`temp'", replace
	restore
	merge 1:1 ticker y  using `temp', update
	drop _merge 
	}

	
* Add risk-free rate to compute excess returns for alternative asset classes	
	preserve
	keep if ticker == "UKbills" | ticker == "UKlongbond"| ticker == "USbills"| ticker == "USTreasury"
	keep y ryreal ticker 
	reshape wide ryreal , j(ticker ) i(y) string
	tempfile temp
	save "`temp'", replace
	restore

	merge m:1 y  using `temp'
	drop _merge 


* restrict sample 
local restrict  `" "drop if inrange(y,1974,1994)" "drop if y<1994" "'

foreach lim of local restrict {
preserve

`lim'

* Compute US-Bonds/Uk bonds
	replace excessBonds = ryreal - ryrealUSTreasury if  ticker == "UScorpbonds"| ticker == "SP"
	replace excessBonds = ryreal - ryrealUKlongbond if  ticker == "FTSE"
	replace excessBonds = ryreal - ryrealUSTreasury if  ticker == "DomEquity"
	replace excessBonds = ryreal - ryrealUSTreasury if  ticker == "DomBonds"
	
	
* Compute excess returns US-Bills
	replace excessBills = ryreal - ryrealUSbills if ticker == "UScorpbonds"| ticker == "SP" | ticker == "USTreasury"
	replace excessBills = ryreal - ryrealUKbills if  ticker == "FTSE"  | ticker == "UKlongbond"
	replace excessBills = ryreal - ryrealUSbills if  ticker == "DomEquity"
	replace excessBills = ryreal - ryrealUSbills if  ticker == "DomBonds"
	replace excessBills = ryreal - ryrealUSbills if  ticker == "safe" & y>=1919
	replace excessBills = ryreal - ryrealUKbills if  ticker == "safe" & y<1919

sort ticker	
	
* Compute excess bond sharpe ratio
	by ticker: egen sd_sharpebonds = sd(excessBonds)
	by ticker: egen excessbonds_average = mean(excessBonds)
	gen sharpebonds = excessbonds_average / sd_sharpebonds
	drop sd_sharpebonds
		
* Compute excess bills sharpe ratio
	by ticker: egen sd_sharpebills = sd(excessBills)
	by ticker: egen excessbills_average = mean(excessBills)
	gen sharpebills = excessbills_average / sd_sharpebills
	drop sd_sharpebills
	
* Compute geometric mean real 
	by ticker: egen geommR = sum( ln( ryreal + 1) )
	by ticker: egen numyears = count(ryreal)
	by ticker: replace geommR = geommR/numyears
	by ticker: replace geommR = exp( geommR ) - 1
	drop numyears

* Compute geometric mean nominal 	
	by ticker: egen geommN = sum( ln( ry + 1) )
	by ticker: egen numyears = count(ry)
	by ticker: replace geommN = geommN/numyears
	by ticker: replace geommN = exp( geommN ) - 1
	drop numyears
	
* Compute SD
	by ticker: egen sd = sd(ryreal)
	
* Compute adjusted sharpe ratio
	gen adj_sharpeBILLS = .
	gen adj_sharpeBONDS = .
	gen Z =.
	levelsof ticker, local(c)
	foreach tick of local c{
	sum ryreal if ticker == "`tick'", det
	local z -1.96 // 95% confbond_IDence interval
	replace Z = (`z' + 1/6 * (`z'^2 - 1)*`r(skewness)' + 1/24 * (`z'^3 - 3*`z')*`r(kurtosis)'  - 1/36*(2*`z'^3 - 5*`z')*`r(skewness)'^2 ) if ticker == "`tick'"
	replace adj_sharpeBILLS = excessbills_average / ( `r(mean)' - Z * `r(Var)' ) if ticker == "`tick'"
	replace adj_sharpeBONDS = excessbonds_average / ( `r(mean)' - Z * `r(Var)' ) if ticker == "`tick'"
	}
	
	
* Compute t-test 
	if "`lim'" == "drop if inrange(y,1974,1994)" {
	mat MTp= J(9,1,.)
	local r = 1
	gen asset = 1 if ticker == "ExtBonds" & (y<1974| y>1994)
	foreach ticker in SP FTSE UScorpbonds DomBonds DomEquity USTreasury UKlongbond USbills UKbills{
	replace asset = 2 if ticker == "`ticker'" & (y<1974| y>1994)
	ttest ryreal, by(asset)
	mat MTp[`r',1] = (real(string(r(p),"%4.2f")))
	local r = 1+ `r'
	}
	}
	 
	if "`lim'" == "drop if y<1994" {
	mat MT2p= J(9,1,.)
	gen asset = 1 if ticker == "ExtBonds" &  y>1994
	local r = 1
	foreach ticker in SP FTSE UScorpbonds DomBonds DomEquity USTreasury UKlongbond USbills UKbills{
	replace asset = 2 if ticker == "`ticker'" & y>1994
	ttest ryreal, by(asset)
	mat MT2p[`r',1] = (real(string(r(p),"%4.2f")))
	local r = 1+ `r'
	}
	}
	
	collapse ryreal ry geomm* sd sharpe* excess* adj*, by(ticker)
	
* -----------------------------------------------------------------------
* Panel A 
	
if "`lim'" == "drop if inrange(y,1974,1994)" {

	* label asset classes
		gen bond_ID = 1 		if 	ticker == "ExtBonds"
		replace bond_ID = 2 	if 	ticker == "SP"
		replace bond_ID = 3 	if 	ticker == "FTSE"
		replace bond_ID = 4 	if 	ticker == "UScorpbonds"	
		replace bond_ID = 5 	if 	ticker == "DomBonds"
		replace bond_ID = 6 	if 	ticker == "DomEquity"	
		replace bond_ID = 7 	if 	ticker == "USTreasury"		
		replace bond_ID = 8 	if 	ticker == "UKlongbond"		
		replace bond_ID = 9 	if 	ticker == "USbills"	
		replace bond_ID = 10 	if 	ticker == "UKbills"	
		replace bond_ID = 11	if 	ticker == "safe"	

	gen t =.
	mat MT= J(15,12,.)
		
		sort bond_ID 
		local r = 1
		
	foreach i in 1 2 3 4 5 6 7 8 9 10 11{

		mat MT[`r',1] = (real(string(100*ryreal[`i'],"%4.2f")))
		mat MT[`r',2] = (real(string(100*ry[`i'],"%4.2f")))
		mat MT[`r',3] = (real(string(100*geommR[`i'] ,"%4.2f")))
		mat MT[`r',4] = (real(string(100*sd[`i'] ,"%4.2f")))
		
		mat MT[`r',6] = (real(string(100*excessBills[`i'],"%4.2f")))
		mat MT[`r',7] = (real(string(sharpebills[`i'],"%4.2f")))
		mat MT[`r',8] = (real(string(adj_sharpeBILLS[`i'] ,"%4.2f")))
		
		mat MT[`r',10] = (real(string(100* excessBonds[`i'],"%4.2f")))
		mat MT[`r',11] = (real(string(sharpebonds[`i'],"%4.2f")))
		mat MT[`r',12] = (real(string(adj_sharpeBONDS[`i'],"%4.2f")))
		
		if `r' == 3 | `r' == 5 | `r' == 8{
		local r = `r' +2
		}
		else {
		local r = 1+ `r'
		}

		}
			

* Continue to create Figure 11
	clear 
	svmat MT, names(col)
	gen bond_ID = _n


	sort bond_ID
	keep if bond_ID == 1| bond_ID == 2| bond_ID== 3| bond_ID == 5| bond_ID== 7 | bond_ID== 8 | bond_ID== 14
	keep bond_ID c1 c7

	gen ticker = ""
	replace	ticker = "ExtBonds" 	if bond_ID == 1
	replace	ticker = "SP" 			if bond_ID == 2
	replace	ticker = "FTSE"  		if bond_ID == 3
	replace	ticker = "UScorpbonds"	if bond_ID == 5
	replace	ticker = "DomBonds" 	if bond_ID == 7
	replace	ticker = "DomEquity"	if bond_ID == 8
	replace ticker = "safe"	 		if bond_ID == 14
	

	sort c1
	gen k =.
	replace k = 1 	if bond_ID == 5
	replace k = 3 	if bond_ID == 14
	replace k = 5 	if bond_ID == 7
	replace k = 7 	if bond_ID == 3
	replace k = 9 	if bond_ID == 1
	replace k = 11 	if bond_ID == 8
	replace k = 13 	if bond_ID == 2


	sort k

* Figure Panel A 

	#delimit ;	
	twoway
	(bar c1 k if bond_ID == 1 , bcolor(maroon)  yaxis(1) yscale(range(0) axis(1)))
	(bar c1 k if bond_ID == 2 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 3 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 5 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 7 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 8 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 14 ,  bcolor(gs11) yaxis(1))
	(scatter c7 k if bond_ID == 1|bond_ID == 2|bond_ID == 3|bond_ID == 5 
					|bond_ID == 7|bond_ID == 8|bond_ID == 14, yaxis(2) 
			mcolor(gs0) lcolor(gs0) connect(direct) lwidth(medthick) yscale(range(0) axis(2)))
	(pcarrowi 7 2 4.5 3 (3), color(black) )	
	(pcarrowi -2.5 10.5 -1 11 (3), color(gs10) )	
	,
		ylabel( -8 -4 0 4 8,nogrid angle(0) labsize(medlarge)  format(%4.0f))
		ylabel(-0.4 -0.2 0 0.2 0.4,nogrid angle(0) labsize(medlarge)  format(%4.1f) axis(2))
		xlabel(1 `" "Corporate" "bonds (US," "S&P index," "since 1900)" "' 
		5 `" "Domestic" "sov. bonds" "(16 adv." "countries," "since 1870)" "' 
		7 `" "UK equities" "(FTSE)" "' 
		9 `" "Sovereign" "bonds" "(in USD or" "GBP, global" "portfolio)" "' 
		13 `" "Equity" "portfolio" "(16 adv." "countries," "since 1870)" "' 
		11 `" "US equities" "(S&P 500)" "' 
		3 `" `""Risk-free""' "gov. bonds" "(UK/US)" "' 
		, 
		
		nogrid angle(0) labsize(small) format(%4.0f)) 
		xtitle("") 
		ytitle("Return, yearly""average,real, in %", size(medlarge) axis(1))
		ytitle("Sharpe ratio", size(medlarge) axis(2))
		xsize(6)
		graphregion(color(white))
		legend( size(large) symysize(*.3)  symxsize(*.3) region(color(white)) row(2) 
		label(1 "Return (real,yearly average; left axis)") 
		label(7 "Sharpe ratio (excess return/SD; right axis)") 
		order(1 7))
		text(4 9 "Sov." "bonds" , color(white) size(small) lcolor(none) fcolor(none) 
		margin(small) box justification(left))
		text(8 2 "{bf:Sharpe Ratio}" , lcolor(none) color(gs0)  fcolor(none) 
		margin(small) box justification(left))
		text(-3 10 "{bf:Return}" , lcolor(none) color(gs9)  fcolor(none) 
		margin(small) box justification(left))
		
	 ;
	#delimit cr		

	graph export "Figure13_PanelA.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure13_PanelA.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure13_PanelA.tif", as(tif) replace

	restore
 }
	
	

* -----------------------------------------------------------------------
* Panel B
	
else {
	gen bond_ID = 1 if 	ticker == "ExtBonds"
	replace bond_ID = 2 if 	ticker == "SP"
	replace bond_ID = 3 if 	ticker == "FTSE"
	replace bond_ID = 4 if 	ticker == "UScorpbonds"	
	replace bond_ID = 5 if 	ticker == "DomBonds"
	replace bond_ID = 6 if 	ticker == "DomEquity"	
	replace bond_ID = 7 if 	ticker == "USTreasury"		
	replace bond_ID = 8 if 	ticker == "UKlongbond"		
	replace bond_ID = 9 if 	ticker == "USbills"		
	replace bond_ID = 10 if 	ticker == "UKbills"	
	replace bond_ID = 11 if 	ticker == "safe"	

	gen t =.
	mat MT2= J(15,12,.)
			
		
	sort bond_ID 
	local r = 1
	foreach i in 1 2 3 4 5 6 7 8 9 10 11{

		mat MT2[`r',1] = (real(string(100*ryreal[`i'],"%4.2f")))
		mat MT2[`r',2] = (real(string(100*ry[`i'],"%4.2f")))
		mat MT2[`r',3] = (real(string(100*geommR[`i'] ,"%4.2f")))
		mat MT2[`r',4] = (real(string(100*sd[`i'] ,"%4.2f")))
		
		mat MT2[`r',6] = (real(string(100*excessBills[`i'],"%4.2f")))
		mat MT2[`r',7] = (real(string(sharpebills[`i'],"%4.2f")))
		mat MT2[`r',8] = (real(string(adj_sharpeBILLS[`i'] ,"%4.2f")))
		
		mat MT2[`r',10] = (real(string(100* excessBonds[`i'],"%4.2f")))
		mat MT2[`r',11] = (real(string(sharpebonds[`i'],"%4.2f")))
		mat MT2[`r',12] = (real(string(adj_sharpeBONDS[`i'],"%4.2f")))
		
	if `r' == 3 | `r' == 5 | `r' == 8{
	local r = `r' +2
	}
	else {
	local r = 1+ `r'
	}

	}


* Continue to create Figure 11
	
	clear 
	svmat MT2, names(col)
	gen bond_ID = _n


	sort bond_ID
	keep if bond_ID == 1| bond_ID == 2| bond_ID== 3| bond_ID == 5| bond_ID== 7 | bond_ID== 8 | bond_ID== 14
	keep bond_ID c1 c7

	gen ticker = ""
	replace	ticker = "ExtBonds" 	if bond_ID == 1
	replace	ticker = "SP" 		if bond_ID == 2
	replace	ticker = "FTSE"  	if bond_ID == 3
	replace	ticker = "UScorpbonds"	 	if bond_ID == 5
	replace	ticker = "DomBonds" 	if bond_ID == 7
	replace	ticker = "DomEquity"	if bond_ID == 8
	replace ticker = "safe"	 		if bond_ID == 14


	sort c1
	gen k =.
	replace k = 1 	if bond_ID == 5
	replace k = 3 	if bond_ID == 14
	replace k = 5 	if bond_ID == 7
	replace k = 7 	if bond_ID == 3
	replace k = 9 	if bond_ID == 2
	replace k = 11 	if bond_ID == 8
	replace k = 13 	if bond_ID == 1

	sort k

	#delimit ;	
	twoway
	(bar c1 k if bond_ID == 1 , bcolor(maroon)  yaxis(1) yscale(range(0) axis(1)))
	(bar c1 k if bond_ID == 2 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 3 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 5 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 7 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 8 ,  bcolor(gs11) yaxis(1))
	(bar c1 k if bond_ID == 14 ,  bcolor(gs11) yaxis(1))
	(scatter c7 k if bond_ID == 1|bond_ID == 2|bond_ID == 3|bond_ID == 5 |bond_ID == 7|bond_ID == 8|bond_ID == 14, yaxis(2) 
			mcolor(gs0) lcolor(gs0) connect(direct) lwidth(medthick) yscale(range(0) axis(2)))
	(pcarrowi 10 3 7 4 (3), color(black) )	
	(pcarrowi 11.5 9 9.5 10 (3), color(gs10) )	
	,
		ylabel(0 4 8 12,nogrid angle(0) labsize(medlarge)  format(%4.0f))
		ylabel(0 0.2 0.4 0.6 0.8,nogrid angle(0) labsize(medlarge)  format(%4.1f) axis(2))
		xlabel(1 `" "Corporate" "bonds (US," "S&P index)" "' 
		5 `" "Domestic" "sov. bonds" "(16 adv." "countries)" "' 
		7 `" "UK equities" "(FTSE)" "' 
		13 `" "Sovereign" "bonds" "(in USD or" "GBP, global" "portfolio)" "' 
		11 `" "Equity" "portfolio" "(16 adv." "countries)" "' 
		9 `" "US equities" "(S&P 500)" "' 
		3 `" `""Risk-free""' "gov. bonds" "(UK/US)" "' 
		, 
		
		nogrid angle(0) labsize(small) format(%4.0f)) 
		xtitle("") 
		ytitle("Return, yearly""average,real, in %", size(medlarge) axis(1))
		ytitle("Sharpe ratio", size(medlarge) axis(2))
		xsize(6)
		graphregion(color(white))
		legend( size(large) symysize(*.3)  symxsize(*.3) region(color(white)) row(2) 
		label(1 "Return (real,yearly average; left axis)") 
		label(7 "Sharpe ratio (excess return/SD; right axis)") 
		order(1 7))
		text(4 13 "Sov." "bonds" , color(white) size(small) lcolor(none) fcolor(none) 
		margin(small) box justification(left))
		text(11 3 "{bf:Sharpe Ratio}" , lcolor(none) color(gs0)  fcolor(none) 
		margin(small) box justification(left))
		text(12 9 "{bf:Return}" , lcolor(none) color(gs9)  fcolor(none) 
		margin(small) box justification(left))
		
	 ;
	#delimit cr		

	graph export "Figure14_PanelB.png", as(png) replace
	graph export "..\figures_tables_for_publication\Figure14_PanelB.eps", as (eps) replace
	graph export "..\figures_tables_for_publication\Figure14_PanelB.tif", as(tif) replace
	
	restore
	}
}


* Put excel 
* ----------------------------------------------------------------
* Define row names 	
	#delimit;
	local row `" "External sovereign bonds (global portfolio)"
				"US equities (S&P 500)"
				"UK equities (FTSE)" 
				""
				"US corporate bonds (S&P AAA, since 1900)"
				""
				"Domestic sovereign bonds (16 adv. countries, since 1870)"
				"Equity portfolio (16 adv. countries, since 1870)"
				""
				"US Treasuries"
				"UK government bonds"
				"US bills (3-m, since 1835)"
				"UK bills (3-month)"
				""'
	;
	#delimit cr
	
	
* Open output file	
	putexcel set "tables.xlsx", sheet("Table 7")  modify
	
* Export row names 
	foreach n in 5 30{
	local i = `n'
	foreach d of local row{ 							
	putexcel B`i'=("`d'") 
	local i = `i' +1
	}
	}	
	
	matselrc MT MTN , r(1/13) 
	matselrc MT2 MTN2 , r(1/13) 
	
	putexcel (C5)=matrix(MTN) 		
	putexcel (C30)=matrix(MTN2) 	

* Format excel cells
	foreach i in 4 29{	
	local r = `i' -1
	$format1  txtwrap top	
	putexcel C`r':C`i' = "Real total returns (yearly, arithm.)", merge top txtwrap hcenter $format2
	putexcel D`r':D`i' = "Nominal total returns (yearly, arithm.)", merge top  txtwrap hcenter $format2
	putexcel E`r':E`i' = "Real total returns (yearly, geom.)", merge top txtwrap hcenter $format2
	putexcel F`r':F`i' = "SD (nom., arithm. return)", merge top txtwrap hcenter $format2
	putexcel H`r':J`r' = "Benchmark: UK/US bills", merge top txtwrap hcenter border(bottom) $format2
	putexcel L`r':N`r' = "Benchmark: UK/US bonds", merge top txtwrap hcenter border(bottom) $format2
	putexcel H`i' = "Excess return"
	putexcel I`i' = "Sharpe ratio"
	putexcel J`i' = "Adj. Sharpe ratio"
	putexcel L`i' = "Excess return"
	putexcel M`i' = "Sharpe ratio"
	putexcel N`i' = "Adj. Sharpe ratio"
	}

	putexcel B2:N2, border(bottom,double) $format2 overwritefmt
	putexcel B18:N18, border(top,double) $format2 overwritefmt	
	putexcel B4:N4, txtwrap hcenter top border(bottom) $format2 overwritefmt
	
	putexcel B27:N27, border(bottom,double) $format2 overwritefmt
	putexcel B43:N43, border(top,double) $format2 overwritefmt	
	putexcel B29:N29, txtwrap hcenter top border(bottom) $format2 overwritefmt
	
	mata													
	b = xl()
	b.load_book("tables.xlsx")
	b.set_sheet("Table 7")
	b.set_column_width(7,7,1) 		//make  column less wbond_IDe
	b.set_column_width(11,11,1) 		//make  column less wbond_IDe
	b.set_column_width(2,2,30) 	
	b.set_row_height(4,4,60) //make  row hight larger
	b.close_book()
	end
	

* Check for t-test 
	* entire period
	mat list MTp
	* modern period
	mat list MT2p




